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ABSTRACT 


This thesis proposes a least-squares error estimator for line-of-sight, direetion of 
arrival-based loealization and a hybrid souree loealization scheme that addresses 
multipath propagation for non-cooperative sources using random arrays of wireless 
sensors. Taking advantage of the dominant reflections, the proposed solution finds the 
location of a signal source by triangulation using the direction of arrival estimations of 
both the line-of-sight and the reflected components. It uses a space division multiple 
access, spread spectrum-based receiver to generate the direction of arrival estimates. The 
time difference of arrival information is used to discriminate between the line-of-sight 
and the non-line-of-sight signals and to associate the incoming multipath signal with the 
corresponding source and reflector pair. In special cases, the proposed scheme is capable 
of solving the association problem spatially without the need for time difference of 
arrival information. Simulation results are included to demonstrate that the proposed 
scheme provides improved estimates by exploiting the non-line-of-sight information. 
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EXECUTIVE SUMMARY 


A wireless sensor network (WSN) is an autonomous, self-organizing network 
without any pre-established infrastructure or centralized administration. WSNs have been 
used for a wide range of applications where, often, the main goal is to monitor a specified 
phenomenon. One important monitoring task which has recently caught the attention of 
WSN researchers is that of locating a signal source by extracting the information 
contained in the received signal. Two primary approaches to source localization have 
evolved; direction of arrival (DOA) based techniques and time difference of arrival 
(TDOA) based techniques 

This thesis proposes a least squares error estimator for DOA localization which is 
unbiased when the noise is Gaussian-distributed with zero mean. This estimator solves an 
over determined Vandermonde system of equations which is known to be 
computationally efficient and accurate. 

Based on this least squares error estimator, this thesis proposes a passive source 
localization scheme which exploits the non-line-of-sight (NLOS) signals from non- 
cooperative sources. The proposed solution is a hybrid DOA/TDOA source localization 
scheme and is comprised of three parts: a DOA estimator, an association algorithm for 
the identified signal bearings, and the source localization scheme itself. The recently 
proposed Space Division Multiple Access (SDMA)-based receiver is used for DOA 
estimation. TDOA information is used to discriminate between the line-of-sight (LOS) 
and the NLOS signals and to associate the incoming multipath signal with the 
corresponding source and reflector pair. In special cases, the proposed scheme is capable 
of solving the association problem spatially without the need for TDOA information. A 
technique is also provided to estimate the position and the orientation of the reflectors 
when site-specific database information is not available. Both centralized and distributed 
variants of the proposed scheme are presented with the latter being of particular interest 
in WSNs. 



The proposed loealization scheme allows a wireless sensor network to (1) perform 
single array localization, (2) perform the localization in a distributed fashion, (3) obtain 
source location estimates with NLOS signals only and (4) improve the location estimates 
compared with those obtained using the LOS information only. 

A simulation model was developed and implemented in MATLAB to evaluate the 
performance of the proposed localization scheme. The simulation results demonstrate that 
the proposed localization scheme provides high accuracy estimates and outperforms the 
LOS-only based localization schemes. This is primarily because more bearings are 
available and the conditioning of the least squares problem is better for the proposed 
scheme. Furthermore, the simulation results also show that the proposed scheme is able 
to provide accurate NLOS-only source location estimates when the LOS paths are not 
available. 
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I. 


INTRODUCTION 


A. INTRODUCTION TO PASSIVE SOURCE LOCALIZATION USING 

WIRELESS SENSOR NETWORKS 

A wireless sensor network (WSN) is an autonomous, self-organizing network 
without any pre-established infrastrueture or eentralized administration [1], WSNs have 
been used for a wide range of applications where, often, the main goal is to monitor a 
specified phenomenon [2], Wireless sensor networks offer numerous advantages when 
compared to traditional wired or wireless networks [3]. Of particular interest, WSNs 
provide greater redundancy since the malfunction of a number of sensors has less impact 
on the overall system performance. Additionally, WSNs can be deployed quickly at low 
cost and are well-suited to use in mobile platforms. Not surprisingly, they have found 
wide-spread interest in emergency and military applications. 

One important monitoring task which has recently caught the attention of WSN 
researchers is that of locating a signal source by extracting the information contained in 
the received signal. That source could be an enemy’s radio location, as in military 
applications, or the direction of arrival (DOA) and the location of a sensed phenomenon 
in an emergency situation like the seismic waves which follow an earthquake [3]. A 
WSN performs the source localization by coordinating the effort of individual sensors 
which act as antenna elements. These sensors are clustered together to form antenna 
arrays which fuse the data collected by the sensors, to carry out the localization task [4]. 

B. RELATED WORK IN PASSIVE SOURCE LOCALIZATION WITH 

WIRELESS SENSORS NETWORKS 

Passive source localization using wireless sensor arrays is a problem addressed 

extensively in existing literature. Two primary approaches to source localization have 

evolved [5]; direction of arrival (DOA) based techniques and time difference of arrival 

(TDOA) based techniques. DOA-based localization systems utilize antenna arrays which 

examine the spatial characteristics of the incident signals to obtain bearing estimates [11], 

[12], [13], [14], [15] and [16]. The bearing estimates are used for the position location 

determination by triangulation [5] [17], [18] and [19]. TDOA localization systems 

1 



estimate the souree loeation using the interseetion of hyperboloids, whieh are the set of 
range differenee measurements between three or more reeeiving sensors. These are 
determined by measurement of the TDOA of a signal between the sensors [2], [6], [7], 
[8], [9], and [10], 

The need for aeeurate loealization requires knowledge of the spatial 
eharaeteristies of the wireless ehannel sinee those eharaeteristies signifieantly affeet the 
performanee of the arrays. Thus, a signifieant ehallenge is the development of realistie 
ehannel models whieh ean aeeurately prediet the behavior of the wireless medium [30]. 
In a real-world deployment, eomplex propagation phenomena lead to uneertainty in 
deeiding whether a direetion of arrival (DOA) eorresponds to the line-of-sight (LOS) 
signal or its refleetion [2]. That uneertainty ean lead to signifieant errors when estimating 
the position of the desired souree. The problem is more severe when a LOS signal does 
not exist as is often the ease in urban environments. Additionally, it has been reported 
that a number of strong refleetions ean be expeeted even in rural areas [31]. These non- 
line-of-sight (NLOS) signals potentially provide additional information that ean be 
exploited in the souree loeation. 

There have been several proposals in literature that eonsider the presenee of 
NLOS signals. The first eategory is eomprised of sehemes whieh attempt to mitigate the 
effeets of the NLOS signals. In [22], the measurements are weighted to emphasize the 
LOS signals, while [23] identifies the arrays that do not reeeive LOS signals and exeludes 
them from the loealization proeess. In both approaehes, the authors seek to minimize the 
impaet of the NLOS signals rather than taking advantage of them. Reeently, a seeond 
eategory of solutions is beginning to emerge that attempts to exploit these NLOS signals. 
The authors of [24] propose a hybrid DOA/TDOA seheme whieh exploits the NLOS 
when the desired souree is eooperative. 

In eontrast, this thesis proposes a passive souree loealization seheme whieh 
exploits the NLOS signals from non-eooperative sourees. It is a DOA/TDOA-based 
seheme whieh uses bearing estimation for loealization through triangulation. The TDOA 
information is used to diseriminate between the LOS and the NLOS signals and to 
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associate the ineoming multipath signal with the eorresponding souree and reflector pair. 
Furthermore, using the NLOS information, the proposed scheme is eapable of performing 
single array loealization. 

C. THESIS OBJECTIVE 

The objective of this researeh is to develop a souree loealization scheme that is 
eapable of non-eooperative souree loealization within the eonstraints of a WSN. The 
existing work on eooperative souree loealization ean then be viewed as a special case of 
this more general solution. 

There are two signifieant eontributions in this work. The first is a proposed least 
squares estimator for DOA-based loealization. The seeond is a passive souree loealization 
seheme whieh exploits the NLOS signals from non-eooperative sourees. The latter is a 
DOA/TDOA-based seheme whieh uses bearing estimation for loealization through 
triangulation. The TDOA information is used to diseriminate between the LOS and the 
NLOS signals and to associate the ineoming multipath signal with the eorresponding 
souree and refleetor pair. Furthermore, using the NLOS information, the proposed 
seheme is eapable of performing single array loealization or NLOS only based 
loealization. Both eentralized and distributed variants of the proposed seheme are 
presented with the latter being of partieular interest in WSNs. 

D. THESIS OUTLINE 

Chapter II provides the baekground to support the proposed work. It introduees 
the fundamental eoneepts of antenna arrays ineluding the response of an array with 
randomly distributed elements. An overview of the wireless environment is also ineluded 
to validate the adopted propagation and reeeived signal models. The ehapter then outlines 
several current approaches to DOA estimation. The reeently proposed Spaee Division 
Multiple Aecess (SDMA)-based receiver is presented and eompared to the eonventional 
MUSIC algorithm. Chapter II eloses with a brief discussion of TDOA estimation. 

The signifieant eontributions of this work are presented in Chapter III. Examining 

LOS-only loealization first, a least-squares estimator for DOA loealization is proposed 

and analysis is provided to investigate the biasing, the impaet of errors, and the 
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conditioning of the proposed estimator. A eomparison between the least-squares and the 
total-least-squares solutions is also presented. To take advantage of a long duration 
signal, a sequential least-squares approaeh is also included. Turning our attention to 
NLOS as well as LOS signals, the proposed loealization seheme which exploits the 
NLOS signals is then deseribed. The ineident signal-souree-refleetor assoeiation 
algorithm is outlined and a teehnique is provided to estimate the position and the 
orientation of the refleetor. Chapter III eoneludes with the souree loealization proeedure 
of the proposed seheme. 

In Chapter IV, simulation results are provided to evaluate the performance of the 
proposed seheme and eompare it to existing LOS-only solutions. 

Chapter V summarizes the significant contributions of this thesis and provides 
some ideas for extending this work in the future. 

Finally, the Appendix ineludes the MATLAB eode used in the simulation. 
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II. BACKGROUND ON SOURCE LOCALIZATION 


This chapter provides the background to support the proposed solution for passive 
source localization of non-cooperative sources using clusters of random or aperiodic 
sensor arrays which form a wireless sensor network (WSN). In this thesis, we consider a 
hierarchical WSN that consists of two levels. The top level is the network of arrays in 
which each array is viewed as a single node. These nodes perform the DOA estimation 
task and the bearings obtained are transmitted to a fusion center, where the localization 
algorithm is executed. In the proposed distributed variant of our solution, each node also 
performs single-array localization. In this case, the source location estimations are then 
simply averaged by the data fusion center. In the second level of the WSN, each array 
contains a number of sensor elements. DOA estimation is conducted at this sensor layer, 
while the source localization is conducted at the array level. An example WSN is shown 
in Figure 1. 



Fusion Center 


Source 


Figure 1. Example WSN comprised of 3 arrays, a single source of interest and a 

fusion center. 
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This chapter begins by discussing the array response to an incident signal. An 
overview of the propagation environment follows and a propagation prediction model is 
adopted to provide a realistic scenario for performance eomparison of the loealization 
schemes. A link budget analysis is then condueted and the signal received by the array is 
computed. Finally, this chapter closes with a discussion of both DO A estimation and 
TDOA estimation. Two methods of DOA estimation are presented, the Multiple Signal 
Classification (MUSIC) DOA estimation method [26] and the new Space Division 
Multiple Aecess (SDMA) reeeiver [25]. A comparison of the angle resolution and the 
aceuraey of both algorithms is included and the SDMA receiver is ehosen as the DOA 
estimation method of the localization process in this thesis. The TDOA estimates will be 
used to assoeiate the multipath components of the received signal with the eorresponding 
souree and the related reflectors. 


A. ARRAY RESPONSE 

The array steering vector contains the responses of all sensors to a souree with a 
single frequency of unit power [26]. The array response varies as a function of direction 
and a steering vector is associated with eaeh direction of interest. The uniqueness of the 
assoeiation is defined by the array configuration [27]. In this section, we begin with a 
discussion of linear arrays followed by the more general case of two-dimensional 
aperiodie and random arrays. 


1. Uniform Linear Array (ULA) 

The configuration of a uniform linear array (ULA) is shown in Figure 2. The 
souree transmits a narrow band signal s(^t) of frequency / and is assumed to be in the far 

field. The array eonsists of N sensors with uniform inter-element spaeingj . With 
respect to the reference node (sensor 1), sensor 2 experiences a time delay of [28] 

d cosd 


At = ■ 


( 1 ) 


where c is the signal propagation speed. The time delay At corresponds to a phase shift 
of the signal equal to 
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( 2 ) 


ls.y/ = In 


d COS0 

I 


This phase shift is the same for every pair of sensors beeause the inter-element spaeing is 
eonstant (uniform). Assuming identieal sensor elements, the steering vector of this array 
is given by 


a[0)- 1 






(3) 


which can in turn be written as 


a 


(») 


1 


- i—d cosO 

e ^ 


-j(V-l)—cos^? 

e ^ 


(4) 


KO 



Figure 2. Uniform linear array of N sensors with inter-element spacing d. The 

source is located in the far field. 

2. Two-Dimensional Aperiodic and Random Arrays 

The geometry of a two-dimensional aperiodic array is shown in Figure 3. The 
sensors of the array are placed in the x-y plane according to some algorithm. Without a 
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loss of generality, the position of the referenee sensor is assumed to be at the origin of the 
coordinate system. The phase difference between sensor i and the reference sensor is 
given as 


Ai// 


In 

T 


[Xj cos0 + y. sm6') 


( 5 ) 


and thus the corresponding steering vector is found to be 


a 


( 0 ) 



{^2 cos^+y2 sin^) 



(x-cos^+y,-sin^) 


( 6 ) 


When the positions of the sensors are chosen by some random process, the aperiodic 
array is known as random array [29]. The steering vector of the random array is identical 


to that of the aperiodic array except that the vector 


a 


(0)= 


1 e 


2 

-}—^{x2 cos^+>'2 sin^) 


is a random vector as in 


7.71 

-y—(x- COS0+5',* sin^) 


( 7 ) 



Figure 3. Two-dimensional array geometry. 


This thesis considers aperiodic, random arrays. Aperiodic sensor arrays have 
several advantages when compared to conventional periodic sensor arrays. Due to the 
non-periodic nature of the sensor spacing, they do not suffer from grating lobes in their 
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spectrum and are not limited to a maximum sensor separation of 0.5 X (0.5 m at 300 
MHz) [29]. This reduces the cost of the construction of such an array, since a smaller 
number of sensors are required. Larger element separation also provides more resilience 
against mutual coupling [29] which occurs between the sensors when one is in the 
vicinity of the other [26]. This coupling effect degrades array performance and is largely 
ignored in most array signal processing literature. Finally, random arrays provide 
flexibility in deployment and can accommodate arbitrary topologies, which are common 
in mobile platforms and WSNs. 

B. SOURCE TO ARRAY 

In this section, we discuss and adopt a radio propagation model of the wireless 
communications channel that will be used in the subsequent performance. Free space 
path loss causes signal strength attenuation of a LOS electromagnetic wave, while with 
NLOS, multipath components can also be attenuated due to reflection, diffraction, and 
scattering [31]. Electromagnetic signals experience attenuation while they travel in space. 
This is the result of spherical energy spread in space. Reflection occurs when a signal 
strikes a surface and is then reflected towards the receiver. Diffraction is the phenomenon 
that occurs when the electromagnetic signal strikes the edge of the corner of a large 
structure compared with the wavelength. Scattering occurs when a signal strikes an object 
that is much smaller compared with the wavelength [32]. These effects lead to complex 
multipath propagation scenarios, especially in areas with a high density of potential 
reflectors and scatterers (e.g., in urban areas). 

1. Propagation Environment 

In this section, we look at signal attenuation in both line-of-sight and reflected 
signals. We present loss expressions for both. 

a. Free Space Loss 

In wireless communications, as the signal propagates through the medium, 
it disperses with distance [33]. This type of attenuation, known as free space loss (L^) , 

can be expressed as the ratio of the signal power between the transmitter and the receiver 
in dB as 
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p ( ^ 

4=101og^ = 201og — 

1 . v /L 


-201og(/l) + 201og(J) + 21.98 dB 


( 8 ) 


where P, is the transmitted signal power, is the received signal power, and d is the 

distance between transmitter and receiver. The free space loss is proportional to the 
square of the distance between the transmitter and the receiver. Thus, as this distance is 
increased, the free space attenuation becomes very large. 


b. Reflection 

When a signal propagating in one medium encounters the boundary of 
another medium, it is partially reflected back to the first medium and partially refracted to 
the new medium [21]. The propagation characteristics of the resulting waves are 
governed by the boundary conditions of the interface. 

A schematic representation of the reflection on a smooth surface is shown 
in Figure 4. The relationships for the angles shown in the figure are given by Snell’s law 
as 

0, = 0„ and 

r- (9) 

sin 6j^ = sin Oj. 


where is the dielectric constant of the material. The attenuation of the reflected signal 
is given by the square of the reflection coefficient as in = Pj |r| . The reflection 

coefficient, T, is a function of the reflected and refracted angles and has a range between 

0 and 1. This coefficient depends on the type of polarization. Thus, for transverse- 
electric (TE) polarized plane, it can be shown to be [21] 

cos cos 0r 

rn = / V —f, (10) 

cos 0j^ + COS 0J, 

while for transverse- magnetic (TM) polarized plane it is [21] 

.cos 0„ - cos Or 
^ 

yjs^ COS 6*^ + cos 0j 
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For s^=5. Figure 5. plots the absolute value of the reflection coefficient 
for both types of polarization as a function of the incident angle. For TM polarization, 
there exists a single incident angle at which Fjj = 0 and no reflection occurs. This angle 
is called the Brewster angle and is expressed as [31] 


sin 


K) = 




( 12 ) 


On the other hand, when the incident angle is 90°, [F^I = |r^ | = 1 and the 

signal is totally reflected for both types of polarization. The path loss because of the 
refection is 


L,=-201og(|r|) 


(13) 



Figure 4. Geometrical representation of the reflection. The incident signal with 
angle 6*, is partially reflected with angle 6*^ and partially refracted with 

angle 6^. 
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Figure 5. |r| as a function of the incident angle for = 5 . 

2 . Radio Propagation Prediction 

Shifting our attention to the multipath effects, the classic radio propagation 
models provide information about the signal power and Doppler shifts of the received 
signal [30]. Flowever, in source localization schemes, the time delay spread and the 
angle-of-arrival spread are also of major importance. A fundamental challenge in 
localization applications is therefore the development of realistic channel models that can 
accurately predict the characteristics of the multipath propagation. These propagation 
models are highly dependent upon the propagation environment and there does not exist a 
single model that covers every environment. Rather, several empirical models have been 
proposed based on measurements from different areas within the specific environment of 
interest. In general, the characterization of these environments is based upon the 
population density and the building architecture (e.g., urban, suburban, and rural). The 
most challenging propagation environment, particularly for source localization problems, 
is the urban environment where the received signal may be dominated by strong 
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reflections and the LOS path is not always present. A typical propagation scenario can be 
seen in Figure 6 and this thesis primarily focuses on this source localization problem in 
urban environments. 



Figure 6. Typical propagation scenario. A LOS path and several NLOS paths are 

shown. 

a. Site-Specific Propagation Prediction 

Most propagation prediction models for urban areas examine the statistical 
parameters of the propagation environment, such as average row spacing and building 
height distributions [21]. That methodology is suitable when the aim of the model is the 
prediction of the average power of the signal received. However, applicability is limited 
in localization applications where the localization scheme is actively attempting to 
exploit the NLOS components of the signal. 

As an alternative, the site-specific prediction model takes advantage of the 
actual mapping of the area under consideration. The mapping can be exploited by using a 
database which contains the footprints of the buildings and any other large objects [21]. 
Those databases must be, in general, three-dimensional. In the special case when the 
assumption of low antennas in tall building environments holds, the databases can be 
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two-dimensional. This is because, in those cases, the primary propagation paths lie 
around the sides of the buildings and no significant propagation paths exist over the 
rooftops [21], 

Site-specific propagation models use ray optical methods and treat each 
reflection as an individual ray in the space. The maximum number of reflections along a 
ray path could be anywhere from six to eight and the number of walls in the database on 
the order of thousands [21] — facts which make the prediction a complex procedure. 
Additionally though (theoretically) the total received power is given by the summation of 
the powers of the individual rays; this strategy does not give a very accurate estimate 
because it ignores the scattered signals. Finally, the accuracy of the databases is on the 
order of 0.5 m in the best case [21], which introduces further uncertainty into the 
prediction. 

All of this notwithstanding, a site-specific approach offers a simplifying 
feature when applied to source localization. When a ray path includes many reflections, 
the power delivered by the ray is reduced dramatically due to the attenuation of each 
reflection and the increased path length. In many cases, only the first reflections are 
strong enough to be exploited. Thus, the power of each individual ray is more important 
in the DOA estimation problem than the total received power and ignoring the 
contribution of the weak components can be an advantage. The localization procedure is 
significantly simplified by ignoring the multi-reflected rays, which could not be exploited 
anyway. 

With respect to the size of the reflector, there is a limitation on the size of 
the region over which reflections occur and ray methods can be used. This can be 
understood by looking at Figure 7 in which a ray from the source is reflected by the 
reflector as shown. The length of the reflector surface is Wg , while 2wp is the width of 
the Fresnel zone, defined by the relative geometry between the receiver and the reflector. 
If Wg cos 6 > 2wp , then the rules for reflection can be used in order to predict the 

propagation characteristics [21]. However, since max{2w^} = VXs where S is the total 
length of the reflected ray in meters, the maximum ray length is then bounded by 
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( 14 ) 


S < —(w^ 008 6^)^ . 

A 

Thus, if ^ = 45° , = 20 m, and A = 0.5 m , (14) limits the ray length to 400 m . If the 

above inequality does not hold, then the object acts as a scatterer [34], 



Figure 7. Fresnel zone and reflector surface characteristics [From 21], 
b. Scattering Prediction 

The propagation rays, as discussed in the previous section, were 
considered discrete lines in space which possess the total power of the corresponding 
multipath component. However, in reality, the power of the signal is distributed around 
these rays which deteriorate the accuracy of the DO A estimation methods. The 
distribution of the power around the ray is the result of the signal scattering caused by 
small objects found around the source. The same phenomenon is observed around the 
reflections which can be considered secondary sources. Thus, the scatterers are grouped 
into clusters, around both the source and the reflectors. 

Several models have been proposed in literature to predict this scattering. 
The distribution of the angles of arrival caused by the scattering can be modeled as 
discrete Gaussian [30] or discrete Laplacian [35] as in 


15 








Gaussian 


( 15 ) 


Pi^[d) = C 2 e Laplacian 

where 6 & \-n,+n^ and the mean value of both funetions is the angle whieh corresponds 

to the associated ray. The parameter a controls the spread of the functions. Typically, a is 
small and the values of both functions are concentrated around the mean. A plot of the 
function for a = 3 is shown in Figure 8. Measurements in [36] indicate that the 
distribution of the incident angle fits the discrete Laplacian function better, since in both 
rural and urban environments, they tend to demonstrate a sharp peak while also 
maintaining long tails. 



Figure 8. Gaussian versus Laplacian scattering functions. 

3. The Adopted Propagation Model 

In summary, the propagation model used in this thesis combines elements from 
the majority of the models discussed above. In a site-specific propagation model, the 
power attenuation of each ray is calculated as the sum of the path loss and the reflection 
attenuation. The site-specific mapping provides the number of reflectors encountered by 
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the ray en route to its destination and the path loss is eomputed using the free spaee path 
loss model for eaeh ray individually. The total power of the reeeived signal is the 
summation of the powers of the individual rays. 

C. RECEIVED SIGNAL 

The reeeived signal model assumes multiple, uneorrelated sourees transmitting 
signals to an N-sensor array. In this seetion, a link budget analysis is performed and the 
reeeived signal-to-noise ratio (SNR) of the ineoming signals is eomputed. Additionally, a 
matrix expression is defined for the reeeived signal as a funetion of the response of the 
reeeiving array. 

I. Link Budget Analysis 

The power level of the ineident signals is eomputed by subtraeting the losses due 
to attenuation, as deseribed previously, from the transmitted power and by adding the 
gain of the array. The gain of the array is given in dB as [29] 

G„=101og(/V) + G,, (16) 

where is the gain of eaeh array element in dB. Expression (16) holds for both ULA 

and aperiodie arrays, although in the latter case, it must be considered an approximation 
[37]. The link budget equation in dB for each signal path is then expressed as 

SNR,, = SNR,,-L,-L,+G„+G, (17) 

where SNR,, is the signal-to-noise ratio available at the receiving array for a specific 
signal path, SNR^^ is the signal-to-noise ratio of the signal transmitted from the 
corresponding source is the free space loss, L, is the path loss because of the 
reflection, G, is the transmitter gain (assumed to be 0 dB in this work) and G,, is the 
array gain. The SNR for signal power P is defined in dB as 
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SNR = 101og - 

^[kTB 


( 18 ) 


where ^is the Boltzmann eonstant equal to 1.38x10^^^ kg K\ T is the system 
noise temperature in degrees Kelvin whieh ineludes both the antenna and reeeiver noise 
and B is the effeetive noise bandwidth of the reeeiver. 

2. Received Signal Model 

Although the signal of eaeh source is considered narrowband in the previous 
discussion, the results can be extended to wideband signals, given the assumption that the 
frequency response of the array is flat over the signal’s bandwidth and the propagation 
time across the array is small when compared to the inverse of the bandwidth [26]. If the 
number of sources is K and each of them transmits a signal 5^ (t), the received signal at 
time t in the array can be expressed as 

x{t) = Y,a{e,)s,{t) + n{t) (19) 

k=\ 

where 

x(t) = [xi(t) x^[t) ... eC^, (20) 

n{t) is the noise vector and a (6(.) is the steering vector for signal s^. (t) with DOA . 
In matrix notation, (20) can be written as 

x{t) = A{^0^~s{t) + n{t) (21) 

where is the array response matrix which contains the responses of each 

sensor for each incident signal and is given as 

A(0] = [a(9,) a(0,) ... a(0,)] (22) 

and ^ is a vector which contains the DOA for each incoming signal. Let L be the 
number of observations with L> K . The received signal can then be written as 
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X =A{@)S + N, 


( 23 ) 


where N e is the noise matrix and S e is the transmitted signal matrix . 

D. DOA ESTIMATION 

This section discusses and compares two methods of DOA estimation, the 
Multiple Signal Classification (MUSIC) estimation and the Space Division Multiple 
Access (SDMA) receiver. 

I. The MUSIC Algorithm 

Most DOA algorithms are based on the computation of the signal correlation 
matrix [32], The received signal correlation matrix and the desired signal correlation 

matrix are defined as 


R:cx=E{x{t)x^ {t)] 

(24) 

II 

(25) 


where H denotes the Hermitian transpose of the matrix and E[a] is the expectation of 

the argument a [26]. If the statistics of the signal and the noise are not known but the 
corresponding processes are ergodic, then the correlation matrices can be approximated 
by averaging a finite number of data observations as [26] 

( 26 ) 

^ l=l 

( 27 ) 

^ l=l 

Both R^^ and R^ are NxN matrices. If additive white Gaussian noise is assumed, the two 
matrices are related as [26] 

C = '4(e)4'»"(e)+CTV (28) 

where is the noise variance and / is an NxN identity matrix. If the incident signals 

are uncorrelated, R^^ is diagonal. When the signals are coherent, is singular [38]. In 
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most cases, the signals are partially correlated and is positive definite. This property is 
very important, sinee, as will be shown later, the DOA estimation algorithms are based 
on the inversion of R^ has N eigenvalues (/{, A 2 , ... and N assoeiated 


eigenvectors £ = [^1 ••• ^n\ which can be obtained by the eigenvalue 

deeomposition. By ordering the eigenvalues from larger to smaller, the eigenvector 
matrix can be divided into two sub-matriees [32] 




E E 


(29) 


These sub-matriees are also called subspaees. E^ has K columns and corresponds to the 

signal subspace, while E^ has N -K columns and corresponds to the noise subspace. An 
alternate way to find the eigenvectors of the autoeorrelation matrix is directly to use the 
reeeived signal matrix X and eigendecompose it by using the singular value 
decomposition (SVD), which gives more stable algorithms [39]. 


The Multiple Signal Classifieation (MUSIC) [26] is the most popular among the 
DOA algorithms based in the subspace decomposition of the eorrelation matrix. The 
desired DO As are estimated by identifying the peaks of the MUSIC spatial 
pseudospectrum, which is given as [40] 




(30) 


The array steering veetors are orthogonal to the noise subspace and, therefore, the peaks 
in the pseudospeetrum represent the DOAs for the desired signal. 

MUSIC can be applied to any arbitrary but known array topology [26] and 
aceordingly, requires accurate array calibration [17]. Furthermore, MUSIC assumes the 
number of sources is known in order to assign the corresponding eigenveetors to the 
signal subspaee and several algorithms have been proposed to do this [41]. Additionally, 
eonventional MUSIC breaks down under near-coherent signal eonditions like those 
which exist when multipath propagation eonditions are present [32]. Again, several 
methods have been proposed to address this, but they are either applicable only to special 
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cases such as spatial smoothing which works for ULA [17] or they are computationally 
intensive such as multidimensional MUSIC [17]. 

2. The SDMA Receiver 

The SDMA receiver is a new method for DO A estimation proposed in [25]. A 
depietion of this reeeiver is shown in Figure 9. The SDMA reeeiver does not rely on the 
subspaee deeomposition of the eorrelation matrix, rather it eross-eorrelates the reeeived 
signal with a pre-eomputed set of array responses for every direetion of interest. 



Figure 9. SDMA reeeiver [From 32]. 


From Figure 9, eonsider an N sensor array. The output of eaeh array element 
u(0 phase-modulated by a set of uneorrelated spreading sequeneesvv^ (t). These 

sequenees ean be any type of pseudorandom [32] or orthogonal [37] sequenee. The 
produeed array outputs are then orthogonal or nearly orthogonal. In matrix notation, the 
output of the array is found to be 

Y = W^X (31) 

where X is given from (22) and W e is the matrix of the spreading sequenees 
whieh is defined as 

IF = [wi(t) w^(t) ... . (32) 
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The signals stored in the virtual array are also modulated by the same set of spreading 
sequences. In matrix notation, the output of the virtual array is given as 

V=W^A^{@), (33) 

where ^ (0) e is the matrix of the array responses of all sensors for all DOAs. The 
correlator cross-correlates the array signal with that of the output of the virtual array as in 

R = V^Y (34) 

where R is a x L matrix. The spatial spectrum of the SDMA receiver is then 

where 

R. = Z^(hO Z^(2,0 ... i^R{K,l) . (36) 

_ /=! /=! /=! 

The peaks of correspond to the DOAs of the incident signals. 

The array used by the SDMA receiver is preferably a two-dimensional random 
array such that the sensor geometry and element phasing is unique for each DOA [32]. 
The spreading technique used further ensures that all received directions of interest are 
uniquely defined. The receiver does not compute the correlation matrix, but just 
correlates the received signal with a pre-computed one. Thus, in contrast to MUSIC, it 
does not rely on the correlation matrix properties. Also, the SDMA receiver does not 
require knowledge of the number of incident signals. Finally, the SDMA receiver does 
not rely on any complex adaptive or slow iterative methods. It “looks” in pre-determined, 
finite set of directions of interest [32] which, by estimating multiple angles 
simultaneously, translates to an instantaneous search through a bank of a finite number of 
expected observations. 

3. Comparison between SDMA Receivers and MUSIC 

A simulation was conducted using MATLAB and a comparison was made of the 
resolution and the accuracy between the SDMA receiver and the MUSIC algorithm. 
Figure 10 shows the correlation magnitude of SDMA and MUSIC for seven incoming 
signals with angular separations of 20 degrees. A random array of 31 sensors, occupying 
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a 25 rectangular area, was used. All the signals had a carrier frequency of 
300 MHz (A = 1 m) with a SNR of 15 dB. The results were averaged across 100 Monte 
Carlo simulations. It is clearly evident that even though the number of signals is large, 
SDMA provides steep lobes in the directions of the incident signals, while MUSIC gives 
main lobes which are difficult to identify. Furthermore, SDMA slightly outperforms 
MUSIC in terms of accuracy as well. Table 1 provides the measured mean value ju and 

the variance <j^ in degrees of the average error of two incident angles. Again the results 
were averaged across 100 Monte Carlo simulations. As expected, both schemes provide 
unbiased estimates since the mean values are near zero, but SDMA has a slightly smaller 
variance. 



Figure 10. Performance comparison of SDMA and MUSIC for random array of 31 
sensors that covers an area of 25 m^. SNR of incident signals is 15 dB. 



A 

(7^ 

MUSIC 

0.0917 

0.6325 

SDMA 

-0.015 

0.5908 


Table 1. Error performance of SDMA and MUSIC. 


23 






























Using the MATLAB histfit.m command, the distribution of the error of both 
schemes was compared with the theoretical zero-mean Gaussian distribution. The results 
are shown in Figure 11 for the SDMA receiver and in Figure 12 for the MUSIC 
algorithm. The results from both schemes can be seen to be closely approximated by the 
error distribution which is indicated by the red curve in both Figures 11 and 12. 


300 



Figure 11. SDMA error distribution. 
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300 



Figure 12. MUSIC error distribution. 


Finally, Figure 13 shows the performanee of the SDMA receiver and MUSIC in 
the case of sparse arrays. The same parameters from Figure 10 were used, except that the 
array area was increased to 2500 m^. A significant improvement in the resolution of both 
methods is observed. As discussed earlier, this is an important advantage when using 
random arrays. The larger spacing between the array elements corresponds to a larger 
array aperture, which provides higher angular resolution [42]. In the sparse-array case, 
SDMA is also seen to outperform MUSIC. 
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Figure 13. Performance comparison of SDMA and MUSIC for a random array of 31 
sensors that covers an area of 2500 m^. SNR of incident signals is 15 dB . 


E. TDOA ESTIMATION 

In TDOA-based localization schemes, the first step estimates the time difference 
of arrival (TDOA) of the incident signal between the sensor nodes in the network. TDOA 
information can be obtained by two general methods: the first involves the subtraction of 
the time of arrival between sensors to obtain the relative difference while the second 
method uses cross-correlation techniques to estimate the desired TDOA [17]. The first 
method requires knowledge of the transmission time which, in the case of non- 
cooperative sources, is not available. Accordingly, only the cross-correlation method will 
be discussed here. 

Assume that the signal, s{t ), is transmitted by an unknown source. Each sensor in 
the cluster will receive amplitude-scaled and time-delayed versions of the transmitted 
signal corrupted by the noise introduced in the channel. Hence, the received signals at 
two different sensors, Vj (?) and .Tj (?) , are given by 
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( 37 ) 


X, (?) = Aj5(?-rj) + nj (?) and 

x,{t) = A,s{t-T,) + n,{t) 

where A, and are the amplitude scaling factors, n,(?) and n 2 (t) are the additive 

noise, and Tj and Tj are the signal offset times at each sensor. Aj and A 2 reside in the 

interval [0,l]. Assumes that the noise is zero-mean Gaussian and the signal and the noise 

are uncorrelated. Equation (37) can be rewritten in the form 

^i(?) = K?) + «i(?) ^ 33 ^ 

^2 (?) = As (? - r) + ^2 (t) 

where, without loss of generality, it is assumed that the first sensor is the one with the 
smaller time of arrival, the amplitudes are normalized by A, and Tj , again without loss of 

generality, is set to 0. The cross-correlation between the signals v,(?) and V 2 (?) can be 
approximated by the estimate 

R^{T) = ^\x,{t)x2{t-T)dt (39) 

where T is the time interval during which the observation is conducted. The lag, t , which 
maximizes (38) is the estimate of the TDOA value [43]. 

F. SUMMARY 

This chapter provided the background to support the proposed solution for passive 
source localization of non-cooperative sources. It began by discussing the array response 
to an incident signal. An overview of the propagation environment followed and a 
propagation prediction model was adopted to provide a realistic scenario for performance 
comparison of the localization schemes. It was shown that through a link budget analysis, 
the signal received by the array can be computed. Finally, the chapter concluded with a 
discussion of both DOA and TDOA estimation techniques. 

In the next chapter we propose a least squares estimator for DOA-based 
localization. Based on this estimator we develop a hybrid DOA/TDOA localization 
scheme which exploits both the LOS and the NLOS signals. 
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III. PASSIVE SOURCE LOCALIZATION USING RANDOM 

SENSOR ARRAYS 

In this chapter, we address the problem of passive source localization and present 
the proposed non-cooperative souree localization scheme. In the context of random 
arrays, we begin with an analysis of LOS-only DOA-based loealization and then move on 
to loealization using both LOS and NLOS signals. 

A. LOS-ONLY DOA-BASED LOCALIZATION 

DOA-based localization schemes use DOA estimates of the incident signals to 
obtain bearings to the position of a source of interest. The bearing estimates are used for 
the position loeation determination by triangulation. In a two-dimensional ease, there are 
two unknowns (the coordinates of the source position). In praetiee, though, more than 
two bearings are required due to finite angular resolution, multipath fading, and noise 
[17]. These bearings form an over-determined system of equations whieh often do not 
interseet at a single point and ean be solved using the least-squares approach. Several 
proposals exist in literature based on the Least Squared Error (LSE) minimization [18], 
[45], while others involve the Minimum Mean Square Error estimation or the Maximum 
Eikelihood [46]. 

I. The Geometry of the Problem 

The proposed least squares estimator for DOA localization is based on the 
geometrie configuration of Eigure 14. Consider K arrays where each array obtains a DOA 
estimate by processing the incoming signal from the source. 6^ is the bearing obtained by 
the array k. We begin by formulating the least squares solution of this problem. 
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Figure 14. LOS DOA-based source localization using three arrays. 


From Figure 14, 


tan6^^ 


y-Jk 
x-x, ’ 


which can be rewritten as 

y-;ctan^, =y,-x,tan^,. 
Doing this for all the arrays yields, in matrix form. 


1 

1 1 

_1 

y 


yj -tan6'jXj 

- tan ^ 2-^2 

1 - tan 9^ 

X 


y,-tan^,x, 


(40) 


(41) 


(42) 


or 


A{9)z=b[9). 


(43) 
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2. Least Squares Solution 

The matrix A(6^) in (42) and (43) has a specific structure and is referred as a 

Vandermonde matrix. The solution of least squares problems involving such a matrix is 
known to be computationally efficient and accurate [47]. 

a. Least-Squares Estimator 

An estimate of the solution to (42) is obtained by minimizing the -norm 
of the residual [48] 

r = b [0) - A(^0 )z . (44) 

This is graphically depicted in Figure 15 where y is the point on the range of A which 
is closest to b . From the normal equations, the solution to (42) can be obtained as 

z=A^{0)b{0) (45) 

where A^(^) = [^A^ (^)A(^)^ is the pseudo-inverse matrix. 



Figure 15. Graphical representation of the least squares problem (From [48]). 
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b. Statistical Analysis of the Least-Squares Estimator 

To show that this estimator is unbiased, let us assume, as discussed 
previously, that the DOA measurements contain a zero-mean Gaussian error SO such that 

o,„ = e+5e. (46) 


Correspondingly, 

\A(>) = A(e)+SA(e) ( 47 ) 


and 

K(s)-b(e)+sb(e). (48) 


Taking the derivative of both and b (6^) from (42), 


dA{e) 

de 


0 

1 



X, 

cos^ 6^ 



cos^ 6*1 

0 

1 


db{6) 

de 

^2 

cos^ 6*2 

and 

cos^ 6*2 

0 

1 




cos^ 6*^ 



cos^ 6*^ 


(49) 


Using the approximations 56 ~ dO, db ~ db and SA^dA, 


0 

1 

1 


56^x^ 

cos^ 6* 


cos^ 6*1 

0 

56, 


56,x, 

cos^ 6*2 

and 5b = 

cos^ 6*2 

0 

56, 


56,X, 

cos^ 6, 


cos^ 6, 


(50) 


Substituting (47) and (48) into (43) 

{A + SAf [b+5b) = {A + 5Af {A + 5A){z+5z) (51) 

where Sz is the error in the source location estimation. Expanding (51) and ignoring the 
second order terms. 
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A^b+A^Sb +SA^b = A^ Az + A^dAz + dA^Az + A^ A5z . 
Rearranging terms, 

A^[db-SAz) + (a^ +SA^)(b- Az) = A^ASz . 

With b - Az = 0, reduees to 

(a^a) ‘ A^ [5b -5Az) = 5z. 

Using the pseudoinverse we finally arrive at 

Sz =A^ [Sb -SAz) 


whieh has the expeetation 

E{Sz} = A^E[Sb-SAz}. 


The veetor inside the expeetation of (53) is 


5b -5Az 


cos^ 9 ^ 
(x-x^) 
cos^ 6^2 


59^ 

59, 


cos^ 9 /^ 


59, 


and has a mean value of 


E[5b-5Az] 


cos^ 9 ^ 
cos^ 9, 


E{59,] 

E{59,] 


[^-^k) 

cos^ 9 ,^ 


E{59,] 


(52) 


(53) 


(54) 


(55) 


However, £'|(5'^^} = 0 foreveryA:e(l, sinee the DO A error is zero- 
mean and, therefore, E[5b -5Az) = 0. From (53), it follows that E[5z] = ^ and the 
proposed estimator is unbiased. 
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The accuracy of the proposed scheme, and of that of any DOA-based 
localization scheme, is affected by the error in DOA estimation. A metric of the accuracy 
of these schemes can be explained in terms of the covariance matrix [49] 

a={A"C,-'Ay (56) 

where Cg e is the covariance matrix of the DOA error estimates (assumed to be 
zero mean white noise sequence in this case). Cg is a diagonal matrix and each element 

along the main diagonal represents the error variance of the measured DOA for the 
corresponding array. If the statistics of the error sequence are known and each array has a 
different noise variance, appropriate weighting of the least-squares solution leads to an 
optimum estimator [49]. The weighting matrix is chosen to be the and the general 
form of the weighted least-squares estimator is 

A..={A'<C,'AyA'‘cAb. (57) 

c. Conditioning of the Least-Squares Solution 

The conditioning of the least-squares problem captures the perturbation 
behavior of the least-squares solution [47]. A least-squares problem is characterized as 
either well-conditioned or ill-conditioned. An ill-conditioned problem is one in which a 
small perturbation in the observed data leads to a large error in the estimated solution and 
the condition number is a metric used to quantify the conditioning. A well-conditioned 
problem has a condition number equal to one. As the condition number increases, the 
problem becomes increasingly ill-conditioned. 

In the least-squares problem of (42), the condition numbers which 
describe the solution z , with respect to the perturbations of the matrix A and the vector 
b as shown in Figure 15, are given by the equations of Table 2 [47] where 

/f(A) = —, l</f(A)<oo , 

6’ = cos ' = , 0<6’<— , 
b 2 
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r| = 


y 


0 <;; < 


and cTj and Oj are the maximum and the minimum singular values of A, respectively. 

Here, the rank of the matrix A is 2 and, therefore, two singular values exist. The 
conditioning of the least-squares problem clearly depends on the geometry of the specific 
scenario under consideration. 



Condition number 

b 

k[A) 

rjcosd 

A 

, . k{A^ iaaO 

/c(A)-t ^ ' - 

11 


Table 2. Condition numbers of the least-squares problem [From 47]. 


d. Geometric Dilution of Precision 

The geometry formed by the arrays which participate in the localization 
procedure and the source of interest affects the accuracy of the estimated source position 
[5]. As the distance between the source and the baseline increases relative to the length of 
the baseline, the accuracy of the least-squares solution decreases. This is because larger 
relative distances correspond to smaller differences in bearings at the arrays. Thus, 
relatively large source distances increase the condition number of matrix A and tend to 
make it singular. This, in turn, affects the stability of the least-squares problem. This 
phenomenon is called geometric dilution of precision and is illustrated in Figure 16. As 
the relative distance of the source from the baseline is increased, the uncertainty area, 
formed by the intersections of the estimated bearings is also increased. Furthermore, the 
angle between the orthogonal bisection of the baseline and the bearing to the source from 
the midpoint of the baseline also affects the accuracy of the solution. As illustrated in 
Figure 17, as angle 0^^ increases, the uncertainty area becomes larger. 


35 




Figure 16. 


Figure 17. 


Array 1 




Array 2 


The effect of distance to the source on the geometric dilution of the 

precision. 


/ 



Array2 


Array 1 



Array2 


The effect of the bearing to the source on geometric dilution of the 

precision. 
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3. 


Total Least Squares 


In the proposed least squares estimator, both A and b can be affected by errors. 
For a least-squares problem, the total least squares (TLS) method is known to compensate 
for errors in matrix A [49]. Let C = [^A with Ce . The singular value 
decomposition of matrix C is 

C = t/Ey" (58) 


where [/e E e and L e . Matrices [/ and V are unitary and E is a 
diagonal matrix of the form 



0 

L, 

0 


0 

0 


(59) 


where its diagonal elements are the singular values of C, satisfying cTj^ > ^ cr^c ■ 

The TLS solution of the least-squares problem using the total least-squares estimator can 
then be shown to be [49] 

^TLS = - C^3cl) ' • (60) 

A simulation was conducted in MATLAB to compare the performance of the 
proposed model using both the LS and the TLS methods. The DOA estimation error was 
assumed to be zero-mean Gaussian with a standard deviation of 0.5 degrees. The results 
are shown in Figure 18 as a function of the number of sensor arrays. The LS solution can 
be seen to outperform the TLS solution for small numbers of arrays, while the TLS 
solution is slightly better for 10 arrays or more. 
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Figure 18. RMS error for LS vs. TLS for 1000 Monte Carlo simulations with 

<^DOA = 0-5° • 

These results are surprising and eontradiet existing literature. As mentioned in 
[49], the TLS is a deregularizing proeedure. The eonditioning of the TLS problem should 
always be worse than the eonditioning of the respeetive LS problem. The ratio of the 
smallest singular value of A to the smallest singular value of C has been defined as a 
metrie to quantify the instability of the TLS and to indieate whether the LS outperforms 

the TLS [49]. As this ratio, in this ease, approaehes unity, the TLS tends to be 


unstable and the LS solution performs better. This ratio was eomputed for the simulation 
of Figure 18 and is shown in Table 3. The results ean be seen to agree with the 
performanee of the TLS algorithm. 


Arrays 

3 

4 

5 

6 

7 

8 

9 

10 

11 

^2 A 

5.4513 

5.4934 

6.3222 

6.9740 

7.S332 

8.9594 

9.9635 

11.0359 

12.0155 












Table 3. Stability measure for the TLS algorithm. 
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4. Sequential Least Squares 


As time progresses, more incoming signal samples may arrive at the array from 
the same source. Hence, more data will be available and, by exploiting it, the estimate of 
the source location can be improved. In this case the localization procedure becomes a 
dynamic task and an algorithm which continuously updates the information matrices is 
needed [48]. If the noise is uncorrelated (which implies that the covariance matrix Cg is 


diagonal), then A can be computed sequentially as 

~A[k-\\ 


A[k] = 


a" [k] 


(61) 


where A[k] is the data matrix of kmeasurements, A[k-l] is the kx2 matrix of the 


previous k-l measurements and a^[k] is the k'*' measurement. Recall that 


—H 

a 


[k] = [l tan6(,] for the proposed scheme. Thus, the sequential estimator is [48] 


^[k] = z[k-l] + D[k](b[k]-a^[k]z[k-l]) 


where 


D[k] = 


Cg [k -l]a [k] 


+ a^ [k][k -1]a [k] 

and the covariance update is given as 

Cg[k] = (l-D[k]a^[k])Cg[k-\]. 


(62) 


(63) 


(64) 


B. LOS AND NLOS DOA-BASED LOCALIZATION 

The proposed hybrid DOA/TDOA LOS and NLOS localization scheme is 
implemented in three steps. A block diagram of the proposed scheme is shown in Figure 
19. The first step implements the DOA estimation using the SDMA receiver as described 
in the previous chapter. A detailed description of steps 2 and 3 (bearing association and 
source localization, respectively) are included in the following subsections. 
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Step 1 


Step 2 


steps 



Figure 19. Block diagram of the proposed LOS and NLOS localization scheme. 

1. Association of Bearings 

Once the DOAs have been determined for the incident signals, they are then 
associated with each other according to source-reflector pairs. Specifically, for each 
source-reflector pair, a LOS and a NLOS is identified based on knowledge of the 
reflector position and geometry. This reflector knowledge can be obtained either from 
static databases which contain the footprint of the large objects within the environment 
under consideration [21] or dynamically by sending a beacon from a known location and 
using the reflector position estimation algorithm which is described in the following 
subsection. When neither of these tools is available, the proposed scheme is capable of 
performing reflector position estimation. In this case, the problem becomes one in which 
the scheme simultaneously performs both reflector mapping and source localization. In 
the remainder of this section, we begin with a discussion of how this reflector mapping is 
accomplished and then present the association scheme. 

a. Reflector Position Estimation 

To discuss the reflector mapping algorithm, consider Figure 20 where a 
single source transmits a narrowband signal and N reflectors generate NLOS signals that 
can be viewed as secondary sources. A minimum of three arrays must be available to 
solve the reflector estimation problem. Let 0. ^ be the bearing associated with array i and 
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reflector n and 6*,^ the LOS bearing from source S to array i. The bearing 


correspondences can be illustrated with a + 4,3A/^ +3)-bipartite graph as shown in 

Figure 21. There are N + 4 vertices in the graph and 3N + 3 edges. The number of edges 
corresponds to the total number of bearings. The interrelationship between the bearings is 
determined in two steps as shown in Figure 22. In the first step, all the intersections 

between the bearings of array 1 and array 2 are found. In this step, a total of (A^ + l)^ 

2x2 systems of equations are solved. In step two, the residuals between those points and 
the bearings of array 3 are computed. The A^ + 1 smaller residuals indicate which 
intersections of array 1 and array 2 bearings correspond to which bearing of array 3. This 

second step involves the solution of (A^ + l)^ equations. Thus, the total number of 


required operations grows as . The source and the image sources related to the 
reflectors are estimated by solving for the interrelated bearings using the least squares 
estimator of the proposed scheme. The orientation and the position of reflector ^ 
according to Figure 20 is then given as 


with 


Zfi 


y + y 


Rn 


2 


(9. = 


0° 



f » ^ 

arctan 





x = X 


Rn 


-90" 


Rn 


(65) 


( 66 ) 
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Figure 20. 


Figure 21. 



Unknown reflector position and orientation estimation with three arrays. 



Representation with a (N + 4,3N + 3) biparte graph of the bearing 
correspondences of the reflector source pairs for N reflectors. 
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Figure 22. Two-step proeedure for finding the bearing eorrespondenees of the 

refleetor souree pairs 


b. Bearing Association using Expected TDOA Estimation 

Onee the positions of the refleetors have been determined, eaeh array ean 
assoeiate its bearings using a set of time markers. From Figure 23 and using the sinusoids 
law, 


‘/wc 




‘LOS 


sin(^^-^J sin(^^ + ^^-26'J sin(2^^-2^J 


(67) 


where is the distanee between the array and its image (with respeet to the refleetor), 
is the orientation of the refleetor, is the LOS bearing and is the NLOS bearing. 
However, 


Ad <7^05 


^INC ^REF 


( 68 ) 


and, eombining (67) and (68), 
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(69) 


AJ 


sin(2^^ -20^)- sin +0^-20^)- sin ^ 
sin(^«+^i- 26 'J 


We also note that 


d 


REF 


d 


AA' 


2sin(^«-6»J 


and 



(70) 


where c is the propagation speed and Ar is the desired expeeted TDOA between the LOS 
and the NLOS signals. Substituting (69) into (70) the expression for Ar is 

(P 


Ar = 


sin (- 2^„) sin ) 


2c 


where 


^ = sin( 26 >^ 




(71) 


(72) 


The TDOAs for all combinations of ordered pairs of bearings are solved 
using (71) for each reflector to give the complete set of possible expected TDOAs which 
can be shown to number 


W 


'm(a + i)^ 




(73) 


TDOAs where M is the number of transmitting sources. It should be emphasized that this 
set is derived using the DOA estimations from the first step of the proposed scheme. 

The set of expected TDOAs is compared with the observed TDOAs, 
which have been calculated from the time markers using a weighted version of the 
generalized cross correlation estimator. For a given source, the bearing association is then 
accomplished by matching ordered pairs of the expected TDOAs to the observed 
TDOAs. 
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Image 



Figure 23. Single-array, single-souree, single-reflector scenario used to illustrate the 
expected TDOA relationship between the LOS and the NLOS signal. 

c. Spatial Bearing Association 

As depicted in Figure 24, when a single source is considered and the exact 
footprint of the reflectors is known, the association problem can be solved spatially 
without the need of TDOA estimation. This saves processing load which is desirable in 
wireless sensor network applications. For the single reflector scenario, two cases can be 
identified. In the first, one bearing is outside of the sector defined between the reflector 
and the array as in Figure 24(a). This bearing can be clearly identified as the LOS signal. 
In the second case, both bearings reside inside the sector as in Figure 24(b). In this case, 
the localization problem (discussed in the next section) is solved twice for each of the 
potential LOS bearings. These solutions will generate two distinct points. One is the 
position of the source while the second is an image of the source related to the specific 
reflector. The latter point can be readily determined to be unrealistic based on the known 
location of the reflector and the former is identified as the actual position of the source. 
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Array 

(a) (b) 

Figure 24. Spatial bearing assoeiation when one source transmits and the exact 

footprint of the reflector is known. 

When multiple reflectors are considered, each one defines an independent 
sector and, in the case where a sector contains two bearings, the association can be solved 
as above. When each sector contains at most one bearing, then the least squares 
localization problem described in the next section is solved as many times as the number 
of the bearings by considering a different bearing as the LOS each time. The least squares 
formulation, which converges to a solution, simultaneously solves both the association 
problem and the localization problem. If no solution results from the procedure, then no 
LOS signal exists and the localization problem is solved by considering all bearings as 
NLOS. 

In the preceding discussion, it is assumed that all reflections are “single 
bounce.” This is a reasonable assumption provided a signal strength threshold is used at 
each receiver. 
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2. Localization using both LOS and NLOS Signals 

The third step in the proposed scheme is the source localization algorithm itself. 
In this section, we begin by describing the single array case and then expand the 
localization procedure to the multiple arrays case. For the latter, we propose both 
centralized and distributed solutions. 


a. Single Array Localization 


As we have noted earlier, the proposed scheme is capable of providing 
single array localization. As shown in Figure 25, if at least one NLOS bearing is 
available, the source position is found as the intersection of the LOS bearing and the 
bearing with angle 29g - 9^ with respect to the image array. This image array is defined 


in a similar manner to the image source described earlier and its position is given by 


1 

_1 


1 

1 

CN 

1_ 

1 

_1 


_2yRA-yA_ 


(74) 


We also note that 


and 


=tan6»„v^ + y^„ 


(75) 


tan (90-6^^) 

Combining (75) and (76), we see that 

1 


— ■ 


1 + tan 0^ 


■Xa+' 


^A-^RA 


tan 6^ 

_ o 

1 + tan^ 9g 


{yA-yRo) 


and 


= 


1 + tan^ 9 ^ 


+ 


tan^ 9^ 

1 + tan^ 9^ 


{yA-yRo)- 


Substituting (77) and (78) into (74), we have 

(l- tan^ 9„)xa + 2 tan6>, () 

_ 1 + tan^ 

2 tan 9^,Xa + (2 tan^ 6>, -1) y^ - 2 tan^ 9j^^ 
1 + tan^ 9^ 


(76) 


(77) 


(78) 


(79) 
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Using (42), the least squares formulation of the single array loealization problem when a 
single NLOS signal exists in eonjunetion with the LOS signal is 


1 -tan 6*^ 

1 -tan(26»„-^^) 


jA-tan^^x^ 

3^A'-tan(2^„-^^)xA, 


(80) 

(81) 


When the refleetions of at least two refleetors are present the loealization 
task ean be performed by using only the NLOS signals. Considering the seenario of 
Figure 25, matrix A and the eorresponding veetor b in (42) beeome 


A = 


1 

1 


-tan( 26 '„i-^«i) 

-tan( 26 '„ 2 -^« 2 ) 


, b = 


yAi.-tan(20„i-^«i)xAi, 

yAi.-tan( 26 '„i-^«i)xAi, 


(82) 


R1 



Figure 25. Single-array NLOS-only loealization with two refleetors. 


b. Proposed Multiple-Array Centralized Localization Scheme 

Consider the multiple array-multiple refleetor seenario of Figure 26 in 
whieh K arrays have available the NLOS refleetions from N refleetors. This problem is 
over-determined and ean be solved using the least-squares estimator proposed earlier. 

The formulation of the matrix ^ and the eorresponding veetor b is 
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(83) 



"A" 




A 


h-y 

A = 


, b= 

2 




■ 1^^ 
_1 


where e the matrix of data obtained by the k"' array and e is the 

corresponding vector. These can be shown to be 



1 - tan 0, 



II 

1 -tan(26»,_-^^J 

and bj, = 

y,, _-tan(26»„__ 


_1 -tan(26»„,^-^,J_ 





Figure 26. Single-source localization using multiple arrays in the presence of 

multiple reflectors. 

This localization scheme is a centralized scheme which involves a data 
fusion center to process the incoming data collected by the arrays. In a WSN, the fusion 
center may be an array of sensors as well. The centralized approach has several 
fundamental drawbacks when implemented in WSNs. The most important is the 
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increased energy eonsumption of the fusion eenter. Thus, the array assigned to aet as the 
fusion eenter will have a shorter lifetime due to the added proeessing and eommunieation 
load whieh will affeet the lifetime of the entire WSN. 

c. Proposed Multiple-Array Distributed Localization Scheme 

As an alternative to the eentralized approaeh above, a distributed solution 
ean be implemented that solves the loealization problem for eaeh single array 
individually and then averages the individual array estimates. A two-step approaeh, this 
solution is summarized in Figure 27. This distributed approaeh is of partieular interest in 
wireless sensor network deployments sinee it distributes the proeessing load among 
multiple sensor nodes and reduees the eommunieation burden whieh is the primary 
energy eonsumer [2]. 



Figure 27. Proposed distributed loealization seheme. 

C. SUMMARY 

In this ehapter, we proposed a least squares estimator for DOA-based loealization 
whieh is unbiased when the noise is Gaussian-distributed with zero mean. This estimator 
solves an over determined Vandermonde system of equations whieh is known to be 
eomputationally effieient and aeeurate. 

Based on this least squares error estimator, we proposed a passive souree 
loealization seheme whieh exploits the NLOS signals from non-eooperative sourees. The 
proposed solution is a hybrid DOA/TDOA souree loealization seheme and is eomprised 
of three parts: a DOA estimator, an assoeiation algorithm for the identified signal 
bearings, and the souree loealization seheme itself The reeently proposed Spaee Division 
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Multiple Access (SDMA)-based receiver was used for DOA estimation. TDOA 
information was used to discriminate between the line-of-sight (LOS) and the NLOS 
signals and to associate the incoming multipath signal with the corresponding source and 
reflector pair. It was shown that, in special cases, the proposed scheme is capable of 
solving the association problem spatially without the need for TDOA information. A 
technique was also provided to estimate the position and the orientation of the reflectors 
when site-specific database information is not available. Both centralized and distributed 
variants of the proposed scheme were presented with the latter being of particular interest 
in WSNs. 

In the next chapter, we evaluate the performance of the proposed passive source 
localization scheme using MATLAB. 
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IV. SIMULATION RESULTS OF THE PROPOSED 
LOCALIZATION SCHEME 

In this chapter, simulation results are provided to validate the performance of the 
proposed localization scheme. We begin with a discussion of the MATLAB simulation 
code, the supporting link budget analysis, and the performance metrics. Following this, 
simulation results are provided and analyzed for both the single reflector and the multiple 
reflector cases. 

A. SET-UP OF THE SIMULATION 

The performance of the proposed localization scheme is validated using the 
MATLAB simulation environment. Details of the underlying simulation are provided in 
the following subsections. 

1. Matlab Simulation 

The MATLAB code used in this thesis is ineluded in the Appendix and is 
implemented in three steps aeeording to the block diagram of Figure 19. The operational 
scenario simulated in MATLAB contains one souree and two reflectors with a 
narrowband signal is transmitted at a given SNR from the source. The propagation 
environment is modeled as diseussed in Chapter II and a link budget analysis is 
conducted. The incident signal at the arrays is processed and the response of each array is 
given by the function resignal.m. The SDMA receiver is simulated with the MATLAB 
function doa sdma.m. Once the DOA estimations have been determined, the simulation 
associates the reflector-source bearing pairs and performs the localization. Finally, the 
MATLAB code simulates and eompares the LOS-only loealization seheme and the 
proposed localization scheme whieh utilizes both LOS and NLOS signals. 

2. Environment Simulation and Link Budget Analysis 

Using the seenario shown in Figure 28, we now provide an example link budget 
analysis. This single-source, single-refleetor example can be expanded to seenarios with 
multiple refleetors or sourees. The souree transmits a narrowband signal with 30 kHz 
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bandwidth at a carrier frequency of 300 MHz . The signal is received at the array via two 
multipath components. The first path is the LOS, while the second is a “single bounce” 
NLOS signal. The power of the transmitted signal is 0.1 mW . The source is located at 
coordinates = (170 m,230 m), while the array is at the origin. The reflector has 

orientation 0^ =10“ with respect to the x-axis. The reference point is located at 
(0m,400m). The angle 6’^=45° corresponds to the DOA of the LOS path, while 
= 79.48° is the DOA of the reflection. The length of the LOS path is 282.84 m , while 
the length of the NLOS path is 667.08 m. Each individual array contains 50 sensors 
randomly distributed in an area of size 15 m x 15 m . 

Assuming a system temperature of 400 K, the transmitted signal SNR is equal 
to 117.81 dB . The received SNRs for the LOS and the NLOS components are 64.93 dB 
and 49.24 dB, respectively. The NLOS path is weaker than the LOS path by 15.69 dB . 
As expected, this is the result of a longer propagation path combined with the reflection 
loss. 



Figure 28. Single-source, single-reflector scenario used in the reported link budget 

analysis. 
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3. 


Performance Metrics 


We model three uncertainties in the simulation scenarios. The first is introduced 
by the noise of the received signal which is considered zero-mean Gaussian. A second is 
generated by the uncertainty in the position of the sensors within the arrays. Finally, the 
different arrays are randomly distributed in the specified area. To capture the 
performance of our proposed scheme, we define a metric based on an estimation 7 of the 
distance of the source from the array r and given as 

/r = (RMS error) 

where 

Y M 

p^ = —^(r-r) (mean error), 

P i=\ 

1 M ^ 

(T^^ = —^(r(error variance) 

P i=i 

and P is the number of Monte Carlo simulations. 

B. SIMULATION RESULTS OF THE PROPOSED LOCALIZATION 
SCHEME 

In this section, simulation results are provided to evaluate the performance of the 
proposed scheme. In the first subsection, we examine the case of a single-source, single- 
reflector. In the second subsection, we extend this to multiple reflectors. 

In the following discussion, the term known reflector implies that the knowledge 
was obtained from existing databases or by a beacon, while the term unknown reflector 
means that the position of the reflector was estimated as described in Chapter III. The 
mapping of the reflector is not free of errors. A standard deviation of 0.25 m was used for 
the position of the known reflector. This value was chosen to agree with the observations 
of[21]. 


(85) 

( 86 ) 
(87) 


I. Single Source-single Reflectors 

In the first scenario, a single reflector provides strong multipath components and a 
single transmitting source with position = (l70 m,230 m) exists in the field. The 
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reflector has orientation 0^ =10° with respect to the x-axis. The reference point of the 
reflector is jr, located at (0,400 m). The WSN covers a square area of 45 mx45 m 

centered at the origin. The general layout of scenario 1 is shown in Figure 29. Seventy 
five Monte Carlo simulations were conducted and a comparison was made between the 
proposed scheme which exploits the NLOS signals and the LOS-only-based scheme. 
Each array was comprised of 50 sensor elements randomly distributed in an area of 
225 m^ 



Figure 29. Scenario 1: Single source located at (170 m,230 m). Single reflector with 
orientation 9^ =10° and reference point at y^j = 400 m . Variable number 
of arrays randomly distributed in a 2025 m^ area. 

The proposed localization scheme outperforms the LOS-only based localization 
as shown in Figure 30. This is because more bearings are available to the proposed 
scheme to formulate the position estimate. Additionally, the condition number of matrix 
A was found to be smaller for the proposed scheme as shown in Table 4. The signals 
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from a common source (i.e., actual source or refleetor) will tend to be elustered in a 
realistie wireless sensor network seenario sinee the relative distanee to the souree of 
interest is typieally large when eompared with the spaeing between neighboring arrays. 
However, the angle separation between those elusters is also relatively large. This is why 
the least squares problem for the proposed seheme is well eonditioned eompared with the 
LOS-only loealization seheme 



Figure 30. Seenario 1; RMS error for both the proposed seheme (LOS and NLOS 
signals) and the LOS-only based loealization seheme in the presenee of 

known refleetors. 


Array 

3 

4 

5 

6 

7 

hA) 

57.8158 

30.3956 

26.7982 

25.3431 

24.4903 

k{A,) 

2.5343 

2.543 

1.7376 

1.6736 

1.7588 


Table 4. Seenario 1: Condition numbers as a funetion of the number of arrays for 

the LOS-only loealization (Aj) and the proposed loealization seheme (A 2 ), 

respeetively. 
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The effect of the conditioning of the least squares problem for both the LOS-only 
and the proposed LOS-NLOS localization schemes is shown in Figure 31 where the RMS 
error is presented as a function of the distance between the source and the WSN. Three 
arrays of 30 sensors each were used. The proposed LOS-NLOS scheme provides accurate 
source position estimates when the distance is large while the LOS-only scheme is 
numerically unstable and inaccurate. We found that in many of the Monte Carlo 
simulations the LS problem of the LOS-only scheme is ill conditioned for source-array 
distances larger than 600 meters. 
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Figure 31. Scenario 1: RMS error as a function of the distance of the source from the 
arrays for both the proposed scheme (LOS and NLOS signals) and the 
LOS-only based localization scheme in the presence of known reflectors 

A comparison between the centralized and the distributed localization scheme is 
shown in Figure 32. Both implementations perform similarly well, but in this particular 
example, the distributed approach is seen to be slightly more accurate. Additionally, for 
the reasons discussed earlier, the distributed approach is the more suitable approach for a 
WSN. 
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Figure 32. Scenario 1: RMS error for both the centralized and the distributed 

configuration of the proposed 


The accuracy of the proposed localization scheme is reduced when the reflector 
position and orientation is unknown. Figure 33. compares the accuracy of the proposed 
scheme for both kn own and unknown reflectors with the LOS-only based localization 
scheme. An iterative approach, which updates the estimated reflector position, improves 
the accuracy of the proposed scheme. Figure 34 presents this improvement for the 
proposed localization scheme as a function of the number of iterations when the 
sequential LS approach is used for five arrays. 
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Figure 33. 


Figure 34. 



Seenario 1; RMS error for the proposed loealization scheme using kn own 
and unknown reflectors and the LOS-only based localization scheme. 



Scenario 1; RMS error as a function of the number of iterations for the 
proposed localization scheme with an unknown reflector using the 
sequential LS approach with five arrays. 
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As discussed earlier, matrix A of the least squares problem includes errors. Thus, 
the proposed localization scheme was also implemented using the TLS solution. This 
TLS solution outperforms the LS solution as shown in Figure 35. This is in contrast with 
what was found for the LOS-only based localization scheme in Chapter III and is the 
result of the well-conditioned nature of the proposed scheme. However, the improved 
location accuracy is achieved at the expense of the larger processing load required by the 
singular value decomposition which makes it difficult to deploy in a WSN. 



Figure 35. 


Scenario 1: RMS error for the proposed localization scheme using least 
squares solution and the total least squares solution. 


2. Single Source-multiple Reflectors 

The second scenario examines the accuracy of the proposed scheme when two 
reflectors are present. The geometric configuration is shown in Figure 36 and is the same 
as in scenario 1, except that an additional reflector exists with orientation 0^ =10° and at 
reference point = -300 m . Again 75 Monte Carlo simulations were conducted and a 

comparison was made between the proposed scheme which exploits the NLOS signals 
and the LOS-only-based scheme. Each array was comprised of 50 sensor elements 
randomly distributed in an area of 225 m^. 
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Figure 36. Scenario 2; Single source located at (170 m,230 m). Two reflectors, both 

with orientation ^ = 10" and reference point 
at = 400 m and = -300 m respectively. A variable number of arrays 

were randomly distributed in a 2025 m^ square area. 

The proposed localization scheme outperforms the LOS-only based localization 
as shown in Figure 37 when the reflectors are known. Furthermore, when analyzing a 
LOS-based scheme, one must consider the case in which the scheme incorrectly identifies 
a NLOS signal to be LOS. As shown in Figure 38, this misidentification causes the LOS- 
only approach to completely break down and return erroneous results. In contrast, the 
proposed scheme, by design, effectively associates the incoming signals. 
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Figure 37. Scenario 2; RMS error for both the proposed scheme (LOS and NLOS 
signals) and the LOS-only based localization scheme in the presence of 
one and two known reflectors 



Figure 38. Scenario 2; RMS error for both the proposed scheme (LOS and NLOS 
signals) and the LOS-only based localization scheme when the latter takes 
one NLOS signal as the LOS signal. 
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The proposed loealization seheme ean also be used for NLOS-only loealization. 
As shown in Figure 39, these estimates are signifieantly improved when more than one 
refleetor is present (as is typieally the ease in a realistie seenario). This is also due to the 
elustered nature of the bearings, as diseussed earlier, whieh affeet the eonditioning of the 
least squares problem. The eondition numbers of the A matrix as a funetion of the 
number of arrays for NLOS-only loealization in the presenee of one and 2 refleetors are 
shown in Table 5. 
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Figure 39. Seenario 2: RMS error for the NLOS-only loealization of the proposed 
loealization seheme in the presenee of one and two known refleetors. 


Array 

3 

4 

5 

6 

7 

k(A,) 

147.7123 

98.8117 

84.1275 

76.5002 

75.3193 

k{A,) 

3.7182 

3.7283 

3.7319 

3.7366 

3.7435 


Table 5. Condition numbers as a funetion of the number of arrays for the NLOS- 

only loealization of the proposed loealization seheme in the presenee of one 
refleetor (Aj) and two refleetors (Aj) , respeetively. 
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C. SUMMARY 

In this chapter, simulation results were provided to validate the performanee of 
the proposed loealization seheme. We began with a diseussion of the MATLAB 
simulation eode and the performanee metries. Following this, simulation results were 
provided and analyzed for both the single refleetor and the multiple refleetor eases. 
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V. CONCLUSION 


This thesis proposed a least squares error estimator. Based on the proposed 
estimator, this thesis developed a passive souree loealization scheme which exploits the 
NLOS signals from non-cooperative sources. The proposed solution is a hybrid 
DOA/TDOA source localization scheme and is comprised of three parts: a DOA 
estimator, an association algorithm for the identified signal bearings, and the source 
localization scheme itself. This scheme allows a wireless sensor network to (1) perform 
single array localization, (2) perform the localization in a distributed fashion, (3) obtain 
source location estimates with NLOS signals only and (4) improve the location estimates 
compared with those obtained using the LOS information only. The proposed solution 
requires the knowledge of the position of potential reflectors which can be obtained 
dynamically using the included reflector mapping algorithm. A reflector knowledge 
updating procedure based on the sequential least squares was used when the reflectors 
were unknown and the accuracy of the proposed scheme was further improved. 
Simulation results were provided to demonstrate the effectiveness of the proposed 
scheme and the conditioning of the least squares problem associated with it and compare 
it to existing solutions. 

A. SIGNIFICANT CONTRIBUTIONS 

There are two significant contributions in this thesis. The first is a least squares 
error estimator for DOA localization which is unbiased when the noise is Gaussian- 
distributed with zero mean. This estimator solves an over determined Vandermonde 
system of equations which is known to be computationally efficient and accurate. 

The second significant contribution of this thesis is a passive source localization 
scheme which exploits the NLOS signals from non-cooperative sources. The simulation 
results showed that the proposed localization scheme provided high accuracy estimates 
and outperformed the LOS-only based localization schemes. This is because more 
bearings are available and the conditioning of the least squares problem is better for the 
proposed scheme 
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B. FUTURE WORK 


Throughout this thesis, it has been assumed that the positions of the sensors are 
known without errors. The presence of errors affects the accuracy of the DO A 
estimations and in turn the accuracy of the proposed localization scheme. A future effort 
may examine methods to reduce the impact of those errors. 

This thesis did not address initial signal of interest (SOI) detection. Since the SOI 
SNR at the receiving sensor network is dependent to a large degree on the characteristics 
of the non-cooperative emitter, future work could include a study of SOI detection, and 
subsequent localization, as a function of received SNR. 

This thesis assumed stationary sources and arrays. This assumption can be 
extended to include platforms that are moving at low speeds. A future effort could also 
examine the performance of the proposed scheme in tracking applications. 

A signal power threshold was used to ensure that all reflections are “single 
bounce.” Future work may solve the association problem of multiple-hop reflections and 
exploit the information contained in those reflections. 

Finally, an analysis which compares the computational cost of both the 
centralized and the distributed configurations of the proposed localization scheme may be 
a subject for a future work. 
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APPENDIX 


This is the Matlab code used in this thesis to evaluate the performance of the 
proposed source localization scheme in the presence of one or two reflectors. 


%THESIS MATLAB CODE 

%Created by Georgios Tsivgoulis 


clear; clc; clf reset; 

%% TRANSMITTED SIGNAL 

f=3E8; 

TB=1/(4*f); 

Nchips=64; 

SPG = 16; 
cp=3E8; 

lambda = cp/f; 
k = 2*pi/lambda; 
tend=SPC*Nchips; 
t=[0:(tend-1)]*TB/(tend-1); 
s=sin(2*pi*f*t); 

S=zeros(1,length(t)) ; 

S(l,:) =S(1,:)+exp(lj*(s)); 

%Signal SNR 
P=10''-4; 

K=1.38*10^-23; 

Tk=400; 

B=3E4; 

SNR=10*logl0(P/(K*Tk*B)) 


% freq of baseband modulation 
% time bandwidth product = 4 
% Number of Chips / Time seq 
% Samples Per Chip 
% speed of light 
% c/f 

% wavenumber 
% time snapshot end 
% discritization of signal 

%incident signal 


% signal power (watt) 
% Boltzman constant 
% system temperature 
% signal bandwidth 
% signal SNR 


%% TARGET POSITION-REFLECTOR SIMULATION 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

0000000000000000000000000000000000000000 

%Simulation of the source and the arrays 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

00000000000000000000000000000000000000000 


targets position 
reflector position 

source distance 
reflector 1 orientation 
reflector 1 reference point 
reflector 2 orientation 
reflector 2 reference point 


tic 

Xt=170; Yt=230; 
stdev=0.25; 
deviation 
T=[Xt Yt]; 

Rt=sqrt (T (1) ''2+T (2) ''2) ; 

thol=10*pi/180; 

YRol=400; 

tho2=-5*pi/180; 

YRo2=-300; 
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c=input( 'give the number of clusters-receivers:' ) 
array_size=input( 'give the size of the array:') 

N=input( 'give the number of sensors/cluster:') 
nruns=input( 'give the number of Monte carlo estimates:') 

% array position 
Xs=zeros(l,c);Ys=zeros(l,c) ; 

%for i=l: (c-1)/2; 

%Ys(2*i)=15*i; 

%Xs(2*i+l)=15*i; 

%end 

for run=l:nruns; % Monte carlo simulation begin from here 

fprintf( 'run:%2.Of\n' ,run); 

YRol=YRol + stdev*randn; 

YRo2=YRo2 + stdev*randn; 
for i=l:c; 

Ys(i)=45*rand; 

Xs(i)=45*rand; 

tan_thL=(Yt-Ys(i))/(Xt-Xs(i)) ; 
thL=atan(tan_thL); 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

oooooooooooooooooooooooooooooo 

%Simulation of the reflector 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

000000000000000000000000000000 


%REFLECT0R1 

XRsl=(1/(1+tan(thol)^2))*Xs(i)+(tan(thol)/(1+tan(thol)^2))*(Ys(i)- 
YRol); 

YRsl=tan(thol)*XRsl+YRol; 

Xsiml=2*XRsl-Xs(i) ; 

Ysiml=2*YRsl-Ys(i) ; 

thAl=2*thol+atan((Ysiml-T(2))/(T(1)-Xsiml)); %DOA REF actual #1 
XR11=50; YRll=tan(thol)*XRll+YRol; 

XR12=200; YR12=tan(thol)*XR12+YRol; 
thll=atan((YRll-Ys(i))/(XRll-Xs(i)))+5*pi/180; 
thl2=atan((YR12-Ys(i))/(XR12-Xs(i)))-3*pi/180; 

%REFLECT0R2 

XRs2=(1/(1+tan(tho2)^2))*Xs(i)+(tan(tho2)/(1+tan(tho2)^2))*(Ys(i)- 
YRo2); 

YRs2=tan(tho2)*XRs2+YRo2; 

Xsim2=2*XRs2-Xs(i) ; 

Ysim2=2*YRs2-Ys(i) ; 

thA2=2*tho2+atan((Ysim2-T(2))/(T(1)-Xsim2)) ; %DOA REF actual #2 
XR21=50; YR21=tan(tho2)*XR21+YRo2; 

XR22=200; YR22=tan(tho2)*XR22+YRo2; 
th21=atan((YR21-Ys(i))/(XR21-Xs(i)))-5*pi/180; 
th22=atan((YR22-Ys(i))/(XR22-Xs(i)))+3*pi/180; 

%% RECEIVED SIGNAL 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

ooooooooooooooooooooooooooo 

%Signal strength simulation 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

ooooooooooooooooooooooooooooo 
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%Array Gain 
GA=10*logl0(N)+2; 


%Assume dipole antenna elements 


%Free space path loss 

dL=sqrt ( (Xt-Xs (1) ) ''2+ (Yt-Ys (1) ) ''2) ; 

LSL=-20*logl0(lambda)+20*logl0(dL)+21.98; 
dAl=sqrt ( (Xt-Xsiml) ''2+ (Yt-Ysiml) ''2) ; 

LSAl=-20*logl0(lambda)+20*logl0(dAl)+21.98; 
dA2=sqrt((Xt-Xsim2)^2+(Yt-Ysim2)^2); 

LSA2=-20*logl0(lambda)+20*logl0(dA2)+21.98; 

%Reflection loss 

thTl=asin(sin(pi/2-thAl)/sqrt(5) ) ; 

Href1=(cos(pi/2-thAl)-sqrt(5)*cos(thTl))/(cos(pi/2- 
thAl)+sqrt(5)*cos(thTl) ) ; 

Lrefl=20*logl0(abs(Hrefl) ) ; 

thT2=asin(sin(pi/2+thA2)/sqrt(5)) ; 

Href2=(cos(pi/2+thA2)- 

sqrt(5)*cos(thT2))/(cos(pi/2+thA2)+sqrt(5)*cos(thT2)) ; 

Lref2=20*logl0(abs(Href2)) ; 

%Received SNR 
SNRL=SNR-LSL+GA; 

SNRA1=SNR-LSA1+Lref1+GA; 

SNRA2=SNR-LSA2+Lref2+GA; 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

ooooooooooooooooooooooooooooooooo 

%Received Signal simulation 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

0000000000000000000000000000000000 

[B1,B2,B3,B] =resignal(thL,thAl,thA2,N,array_size,k); 
XL=B1*S; %Received signal 
XL=awgn(XL,SNRL, 'measured' ); 

XA1=B2*S; %Received signal 
XAl=awgn(XAl,SNRAl, 'measured' ); 

XA2=B3*S; %Received signal 
XA2=awgn(XA2,SNRA2, 'measured' ); 

X=XA1+XA2+XL; 

%% DOA ESTIMATION 
thESTsdma=doa_sdma(X,B,N); 

%thESTsdma2=doa_sdma(XAl, B, N) ; 

%thESTsdma3=doa_sdma(XA2, B, N) ; 

%thESTsdma=[thESTsdmal;thESTsdma2;thESTsdmaS]; 

%% BEARING ASSOCIATION 

%LOS=find(thESTsdma>thll | thESTsdma<thl2 |thESTsdma>th21 | 
thESTsdma>th22 ); 

REFl=find(thESTsdma<thll & thESTsdma>thl2); 

REF2=find(thESTsdma>th21 & thESTsdma<th22); 

LOS=find(thESTsdma<thl2&thESTsdma>th22) ; 
thLest(1)=thESTsdma(LOS); 
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thAestl(i)=thESTsdma(REFl) ; 
thAest2(i)=thESTsdma(REF2); 
end 

%% LOCALIZATION 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

ooooooooooooooooooo 

%Matrix Formulation 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

0000000000000000000 

for j=3:c; 

for i=l:j 

A(i,:) = [l -tan(thLest (i))]; 

A(i+j,:)=[1 -tan(2*thol-thAestl(i))]; 
A(i+2*j,:)=[l -tan(2*tho2-thAest2(i))]; 
b (i,1)=Ys(i)-tan(thLest(i))*Xs(i) ; 
b(i+j,l)= Ysiml-tan(2*thol-thAest1(i))*Xsiml; 
b(i+2*j,l)= Ysim2-tan(2*tho2-thAest2(i))*Xsim2; 

%Matrix formulation for image points 

Arefl(i,:)=[1 -tan(thAestl(i))]; 

bref1(i,1)=Ys(i)-tan(thAestl(i))*Xs(i); 

Aref2(i,:)=[1 -tan(thAest2(i))]; 
bref2 (i,1)=Ys(i)-tan(thAest2(i))*Xs(i) ; 
end 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

000000000000000000000000000000000000 

%l)LOS-only localization 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

00000000000000000000000000000000000 

Xestl =A(1:j, :)\b (1:j,1) 

Rtestl=sqrt(Xestl(1)^2+Xestl(2)^2) ; 

Y1=A(1:j,:)*Xestl; 
kappal=cond(A(1:j,:)); 

thetal=asin(norm(b(l:j,l)-Yl)/norm(b(l:j,l))); 
etal=norm(A(1:j , :))*norm(Xest1)/norm(Y1); 

Condbxl(run,j-2)=kappal/(etal*cos(thetal)); 

CondAxl(run,j-2)=kappal+(kappal^2*tan(thetal))/etal; 
Yt_sdmal(run,j-2)=Xestl(1) ; 

Xt_sdmal(run, j-2)=Xestl(2 ) ; 

%Bias of the Target location 
biasYt_sdmal(run,j-2)=T(2)-Xestl(1); 
biasXt_sdmal(run,j-2)=T(1)-Xestl(2); 
biasRt_sdmal(run,j-2)=Rt-Rtestl; 

%REFLECTOR DETERMINATION% 

%REFLECT0R1 

Xestref1 =Aref1(l:j,:)\brefl(l:j,l); 

tan_thMl=(Xestref1(1)-Xestl(1))/(Xestl(2)-Xestref1(2)) ; 
tho_est1=((pi/2)-atan(tan_thMl)); 
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tho_est1*180/pi; 

YRo_estl=(Xestl(1)iXestref1 (1))/2- 

tan(tho_estl)*((Xestl(2)tXestrefl(2))/2) ; 

Xestref2 =Aref2\bref2; 

%REFLECTOR2 

tan_thM2=(Xestref2(1)-Xestl(1))/(Xestl(2)-Xestref2(2)); 

tho_est2=((pi/2)-atan(tan_thM2)); 

tho_est2*180/pi; 

YRo_est2=(Xestl(1)tXestref2 (1)) 12 - 

tan(tho_est2)*((Xestl(2)+Xestref2(2))/2); 


XRsestl=(1/(1+tan(tho_estl)*2))*Xs(i)+(tan(tho_estl)/(1+tan(tho_estl)^2 
))*(Ys(i)-YRo_est1); 

YRsestl=tan(tho_est1)*XRsestl+YRo_est1; 

Xsimestl=2*XRsestl-Xs(i) ; 

Ysimestl=2*YRsestl-Ys(i) ; 

XRsest2=(1/(1+tan(tho_est2)*2))*Xs(i)+(tan(tho_est2)/(1+tan(tho_est2)^2 
))*(Ys(i)-YRo_est2); 

YRsest2=tan(tho_est2)*XRsest2+YRo_est2; 

Xsimest2=2*XRsest2-Xs(i) ; 

Ysimest2=2*YRsest2-Ys(i) ; 

%Matrix Formulation 
for i=l:j 

Ar(i,:) = [l -tan(thLest (i))]; 

Ar(i+j,:)=[1 -tan(2*tho_est1-thAest1(i))]; 

Ar(i+2*j, :) = [1 -tan(2*tho_est2-thAest2(i))]; 
br(i,1)=Ys(i)-tan(thLest(i))*Xs(i) ; 

br(i +j , 1)= Ysimestl-tan(2*tho_est1-thAest1(i))*Xsimest1; 
br(i+2* j,1)= Ysimest2-tan(2*tho_est2-thAest2(i))*Xsimest2; 
end 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

00000000000 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

ooooooooooo 

%SCENARIO 1 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

000000000000 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

oooooooooooo 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%2)Scenario 1 : Centralized proposed localization ( 1 known reflector) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 9 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

A2=A(l:2*j, :); b2=b(1:2*j, : ) ; 

Xest2=A2\b2 
Y2=A2*Xest2; 

Rtest2=sqrt(Xest2(1)^2+Xest2 (2)^2) ; 

Yt_sdma2(run,j-2)=Xest2(1); 

Xt_sdma2(run,j-2)=Xest2(2) ; 

kappa2=cond(A2); 

theta2=asin(norm(b2-Y2)/norm(b2)) ; 


73 


eta2=norm(A2)*norm(Xest2)/norm(Y2) ; 

Condbx2(run,j-2)=kappa2/(eta2*cos(theta2)) ; 

CondAx2 (run, j-2 ) =kappa2+ (kappa2''2*tan (theta2 ) ) /eta2 ; 

%Bias of the Target location 
biasYt_sdma2(run,j-2)=T(2)-Xest2(1) ; 
biasXt_sdma2(run, j-2)=T(1)-Xest2(2 ) ; 
biasRt_sdma2(run,j-2)=Rt-Rtest2 ; 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%3)Scenario 1 : Distributed propossed localization (1 known reflector) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

for i=l:j 

A3=[A(i, :);A(i +j, : ) ] ; 
b3=[b(i);b(i + j)] ; 

Test3(:,i)=A3\b3; 
end 

Xest3=mean(Test3, 2) ; 

Yt_sdma3(run,j-2)=Xest3(1) ; 

Xt_sdma3(run, j-2)=Xest3(2 ) ; 

Rtest3=sqrt (Xest3 (1) ''2+Xest3 ( 2 )^ 2 ) ; 

%Bias of the Target location 
biasYt_sdma3(run, j-2)=T(2)-Xest3(1) ; 
biasXt_sdma3(run,j-2)=T(1)-Xest3(2) ; 
biasRt_sdma3(run,j-2)=Rt-Rtest3; 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%4) Scenario 1 : Centralized proposed localization (1 unknown reflector) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

A4=Ar(1:2*j, : ) ; 
b4=br(l:2*j, :) ; 

Xest4 =A4\b4; Y4=A4*Xest4; 

Rtest4=sqrt(Xest4(1)^2+Xest4 (2)^2); 

Yt_sdma4(run,j-2)=Xest4(1); 

Xt_sdma4(run,j-2)=Xest4(2) ; 


%Bias of the Target location 
biasYt_sdma4(run,j-2)=T(2)-Xest4(1) ; 
biasXt_sdma4(run,j-2)=T(1)-Xest4(2) ; 
biasRt_sdma4(run,j-2)=Rt-Rtest4; 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%5)Scenario l:NLOS-only localization (1 known reflector) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

A5=A(j + l:2*j, :) ; 
b5=b(j + l:2*j,l) ; 

XestS =A5\b5; 

Rtest5=sqrt (XestS (1) ''2+Xest5 ( 2 )^ 2 ) ; 

Y5=A5*Xest5; 
kappa5=cond(AS); 

thetaS=asin(norm(bS-YS)/norm(bS)) ; 
etaS=norm(AS)*norm(XestS)/norm(YS); 


74 


CondbxS(run,j-2)=kappa5/(eta5*cos(thetaS)); 

CondAxS(run,j-2)=kappa5+(kappa5^2*tan(thetaS))/eta5; 

Yt_sdma5(run,j-2)=Xest5(1); 

Xt_sdma5(run,j-2)=Xest5(2 ) ; 

%Bias of the Target location 
biasYt_sdma5(run, j-2)=T(2)-XestS(1) ; 
biasXt_sdma5(run,j-2)=T(1)-XestS(2 ) ; 
biasRt_sdma5(run, j-2)=Rt-Rtest5; 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%6) Scenario 1 : Centrilized proposed localization (TLS, 1 reflector) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

nl = size(A4,2); 

C = [A4 b4]; 

[U1 Si VI] = svd(C,0); 

Vll = VI(1:nl,1+nl:end); 

V12 = VI(1+nl:end,1+nl:end); 

Test6 = -V11/V12; 

Rtest6=sqrt(TestG(1)^2+Test6(2)^2) ; 

Yt_sdma6(run,j-2)=Test6(l); 

Xt_sdma6(run,j-2)=Test6(2) ; 

%Bias of the Target location 
biasXt_sdma6(run,j-2)=T(2)-TestG(1) ; 
biasYt_sdma6(run, j-2)=T(1)-TestG(2) ; 
biasRt_sdmaG(run,j-2)=Rt-RtestG; 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

000000000000000000000000000000000000000000000000000 

%7 ) Scenario2 : Centralized proposed localization (2 known reflector) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

oooooooooooooooooooooooooooooooooooooooooooooooooo 

Xest7=A\b 

Y7=A*Xest7; 

Rtest7=sqrt (Xest7 (1) ''2+Xest7 (2) ''2) ; 

Yt_sdma7(run,j-2)=Xest7(1) ; 

Xt_sdma7(run,j-2)=Xest7(2); 

kappa7=cond(A); 

theta7=asin(norm(b-Y7)/norm(b)) ; 
eta7=norm(A)*norm(Xest7)/norm(Y7 ) ; 

Condbx7(run,j-2)=kappa7/(eta7*cos(theta7)); 

CondAx7(run,j-2)=kappa7+(kappa7^2*tan(theta7))/eta7; 

%Bias of the Target location 
biasYt_sdma7(run,j-2)=T(2)-Xest7(1) ; 
biasXt_sdma7(run,j-2)=T(1)-Xest7(2) ; 
biasRt_sdma7(run,j-2)=Rt-Rtest7 ; 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%8)Scenario 2 : Centralized proposed localization (2 unknown reflector) 
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9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

XestS =Ar\br; 

Rtest8=sqrt (XestS (1) ''2+Xest8 (2) ''2) ; 

Yt_sdma8(run,j-2)=Xest8(1) ; 

Xt_sdma8(run,j-2)=Xest8(2); 


%Bias of the Target location 
biasYt_sdma8(run, j-2)=T(2)-XestS(1) ; 
biasXt_sdma8(run,j-2)=T(1)-XestS(2) ; 
biasRt_sdma8(run,j-2)=Rt-Rtest8 ; 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%9)Scenario 2:Distributed proposed localization (2 known reflector) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

for i=l:j 

A9=[A(i, :);A(i + j, :);A(i+2*j, :)] ; 
b9=[b(i);b(i + j);b(i+2*j) ] ; 

Test9(:,i)=A9\b9; 
end 

Xest9=mean(Test9,2); 

Yt_sdma9(run,j-2)=Xest9 (1); 

Xt_sdma9(run, j-2)=Xest9(2); 

Rtest9=sqrt(Xest9(1)^2+Xest9(2)^2) ; 

%Bias of the Target location 
biasYt_sdma9(run,j-2)=T(2)-Xest9(1); 
biasXt_sdma9(run,j-2)=T(1)-Xest9(2); 
biasRt_sdma9(run,j-2)=Rt-Rtest9; 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%10)Scenario 2:NL0S-only (2 known reflector) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

A10=A(j + l:3*j, : ); 
blO=b(j+1:3*j,1); 

XestlO =A10\bl0; 

RtestlO=sqrt(XestlO(1)^2+XestlO (2)^2) ; 

Y10=A10*Xestl0; 
kappalO=cond(AlO); 

thetalO=asin(norm(blO-YlO)/norm(blO)) ; 
etalO=norm(AlO)*norm(XestlO)/norm(YIO); 

CondbxlO(run,j-2)=kappalO/(etalO*cos(thetalO)) ; 

CondAxlO (run, j-2) =kappalO + (kappal0''2*tan (thetalO) ) /etalO; 

Yt_sdmalO(run,j-2)=XestlO(1) ; 

Xt_sdmalO(run, j-2)=XestlO(2); 

%Bias of the Target location 
biasYt_sdmalO(run,j-2)=T(2)-XestlO(1); 
biasXt_sdmalO(run,j-2)=T(1)-XestlO(2 ) ; 
biasRt_sdmalO(run, j-2)=Rt-Rtest10 ; 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%ll)LOS-only localization (fault angle) 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 
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for i=l:j-1; 

All(i+1,:)=[1 -tan(thLest(i+1))]; 
bll(i+l,:)=[l -tan(thLest(itl))]; 
end 

All(l,:) = [l -tan(thAestl (1))]; 
bll(l,:) = [l -tan(thAestl (1))]; 

Xestll =All\bll; 

Rtestll=sqrt (Xestll (1) ''2+Xestll (2) ''2) ; 

Yt_sdmall(run,j-2)=Xestll(1) ; 

Xt_sdmall(run, j-2)=Xestll(2) ; 

%Bias of the Target location 
biasYt_sdmall(run, j-2)=T(2)-Xestll(1) ; 
biasXt_sdmal1(run,j-2)=T(1)-Xestll(2); 
biasRt_sdmal1(run,j-2)=Rt-Rtest11; 

end 

end 

%% ERROR ANALYSIS 

%1)Statistics for LOS-only localization 

bias_sdmal=[biasXt_sdmal' biasYt_sdmal' biasRt_sdmal']; 

bias 

meanbiasRt_sdmal=mean(biasRt_sdmal); 
position error 

varbiasRt_sdmal=var(biasRt_sdmal); 
of position error 

rmsbiasRt_sdmal=sqrt(meanbiasRt_sdmal.^2+varbiasRt_sdmal); 
position error 

meanCondbxl=mean(Condbxl) 
meanCondAxl=mean(CondAxl) 

%2)Statistics for Centralized proposed localization (1 known 
bias_sdma2=[biasXt_sdma2' biasYt_sdma2' biasRt_sdma2']; 
bias 

meanbiasRt_sdma2=mean(biasRt_sdma2); 
position error 

varbiasRt_sdma2=var(biasRt_sdma2); 
of position error 

rmsbiasRt_sdma2=sqrt(meanbiasRt_sdma2.^2+varbiasRt_sdma2); 
position error 

meanCondbx2=mean(Condbx2) 
meanCondAx2=mean(CondAx2) 


%3 ) Statistics for Distributed proposed localization ( 1 known 
bias_sdma3=[biasXt_sdma3' biasYt_sdma3' biasRt_sdma3']; 
bias 

meanbiasRt_sdma3=mean(biasRt_sdma3); 
position error 


%Matrix of 
%Mean 
%Variance 
%RMS 

reflector) 
%Matrix of 

%Mean 

%Variance 

%RMS 


reflector) 
%Matrix of 

%Mean 
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%Variance 


varbiasRt_sdma3=var(biasRt_sdma3); 
of position error 
rmsbiasRt_sdma3=sqrt(meanbiasRt_sdma3.''2+varbiasRt_sdma3); %RMS 

position error 

%4)Statistics for Centralized proposed localization (1 unknown 
reflector) 

bias_sdma4=[biasXt_sdma4' biasYt_sdma4' biasRt_sdma4']; %Matrix of 

bias 

meanbiasRt_sdma4=mean(biasRt_sdma4); %Mean 

position error 

varbiasRt_sdma4=var(biasRt_sdma4); %Variance 

of position error 

rmsbiasRt_sdma4=sqrt(meanbiasRt_sdma4. ''2+varbiasRt_sdma4); %RMS 

position error 

%5)Statistics for NLOS-only localization (1 known reflector) 
bias_sdma5=[biasXt_sdma5' biasYt_sdma5' biasRt_sdma5']; 
bias 

meanbiasRt_sdma5=mean(biasRt_sdma5); 
position error 

varbiasRt_sdma5=var(biasRt_sdma5); 
of position error 

rmsbiasRt_sdma5=sqrt (meanbiasRt_sdma5 . ''2+varbiasRt_sdma5) ; 
position error 

meanCondbx5=mean(CondbxS) 
meanCondAx5=mean(CondAxS) 

%6)Statistics for Centrilized proposed localization (TLS, 
bias_sdma6=[biasXt_sdma6' biasYt_sdma6' biasRt_sdma6']; 
bias 

meanbiasRt_sdma6=mean(biasRt_sdma6); 
position error 

varbiasRt_sdma6=var(biasRt_sdma6); 
of position error 

rmsbiasRt_sdma6=sqrt(meanbiasRt_sdma6. ''2+varbiasRt_sdma6) 
position error 

%7)Statistics for Centralized proposed localization(2 known 
bias_sdma7=[biasXt_sdma7' biasYt_sdma7' biasRt_sdma7']; 
bias 

meanbiasRt_sdma7=mean(biasRt_sdma7); 
position error 

varbiasRt_sdma7=var(biasRt_sdma7); 
of position error 

rmsbiasRt_sdma7=sqrt(meanbiasRt_sdma7.^2+varbiasRt_sdma7); 
position error 

meanCondbx7=mean(Condbx7) 
meanCondAx7=mean(CondAx7) 

%8)Statistics for Centralized proposed localization (2 known reflector) 


reflector) 
%Matrix of 

%Mean 

%Variance 

%RMS 


1 reflector) 
%Matrix of 

%Mean 

%Variance 

%RMS 


%Matrix of 
%Mean 
%Variance 
%RMS 
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%Matrix of 


bias_sdma8=[biasXt_sdma8' biasYt_sdma8' biasRt_sdma8']; 
bias 

meanbiasRt_sdma8=mean(biasRt_sdma8); %Mean 

position error 

varbiasRt_sdma8=var(biasRt_sdma8); %Variance 

of position error 

rmsbiasRt_sdma8=sqrt(meanbiasRt_sdma8. ''2+varbiasRt_sdma8); %RMS 

position error 

%9)Statistics for Centralized proposed localization(2 unknown 
reflector) 

bias_sdma9=[biasXt_sdma9' biasYt_sdma9' biasRt_sdma9']; %Matrix of 

bias 

meanbiasRt_sdma9=mean(biasRt_sdma9); %Mean 

position error 

varbiasRt_sdma9=var(biasRt_sdma9); %Variance 

of position error 

rmsbiasRt_sdma9=sqrt(meanbiasRt_sdma9.^2+varbiasRt_sdma9); %RMS 

position error 

%10)Statistics for NLOS-only localization (1 known reflector 
bias_sdmalO=[biasXt_sdmalO' biasYt_sdmalO' biasRt_sdmalO']; 
bias 

meanbiasRt_sdmalO=mean(biasRt_sdmalO); 
position error 

varbiasRt_sdmalO=var(biasRt_sdmalO); 
of position error 

rmsbiasRt_sdmalO=sqrt(meanbiasRt_sdmalO. ''2+varbiasRt_sdmalO) 
position error 

meanCondbxlO=mean(CondbxlO) 
meanCondAxlO=mean(CondAxlO) 

%11)Statistics for LOS-only localization (fault angle) 
bias_sdmall=[biasXt_sdmall' biasYt_sdmal1' biasRt_sdmall']; 
bias 

meanbiasRt_sdmall=mean(biasRt_sdmal1); 
position error 

varbiasRt_sdmall=var(biasRt_sdmal1); 
of position error 

rmsbiasRt_sdmall=sqrt(meanbiasRt_sdmal1.^2+varbiasRt_sdmal1 
position error 


toe 

SN=3:7; 
figure(1) 

plot(SN,rmsbiasRt_sdmal, 'b—^ ,SN,rmsbiasRt_sdma2, 'r-*') 
xlabel( 'Number Arrays' , 'Fontsize' ,12) 
ylabel( '\gamma,RMS(meters)' , 'Fontsize' ,12) 

legend (' LOS-only localizationProposed localization (known 
reflector)' ) 


%Matrix of 
%Mean 
%Variance 
); %RMS 


) 

%Matrix of 
%Mean 
%Variance 
; %RMS 
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figure(2) 

plot(SN,rmsbiasRt_sdma2, 'g—^ ,SN,rmsbiasRt_sdma3, 'r-*') 
xlabel (' Number of ArraysFontsize 12) 
ylabel( '\gamma-RMS(meters)' , 'Fontsize' ,12) 

legend (' Centralized proposed localizationDistributed proposed 
localization ' ) 

figure(3) 

plot(SN,rmsbiasRt_sdmal, 'b-^' ,SN,rmsbiasRt_sdma2, 'g- 

o' ,SN,rmsbiasRt_sdma4, 'r-*' ) 

xlabel (' Number of ArraysFontsize ', 12) 

ylabel( '\gamma,RMS(meters)' , 'Fontsize' ,12) 

legend (' LOS-only localizationProposed localization (known 
reflector)' , ... 

'Proposed localization (unknown reflector) ') 
figure(4) 

plot(SN,rmsbiasRt_sdma4, 'r—^' ,SN,rmsbiasRt_sdma6, 'r-*' ) 
xlabel (' Number of ArraysFontsize ', 12) 
ylabel( '\gamma,RMS(meters)' , 'Fontsize' ,12) 

legend (' LS proposed scheme(unknown reflector) ', 'TLS proposed 
scheme(unknown reflector)' ) 

figure(5) 

plot(SN,rmsbiasRt_sdmal, 'b—^ ,SN,rmsbiasRt_sdma2, 'r- 

*' ,SN,rmsbiasRt_sdma7, 'r-*' ) 

xlabel( 'Number Arrays' , 'Fontsize' ,12) 

ylabel( '\gamma,RMS(meters)' , 'Fontsize' ,12) 

legend (' LOS-only localizationProposed localization (1 known 
reflector)' , ... 

'Proposed localization (2 known reflector) ') 
figure(6) 

plot(SN,rmsbiasRt_sdmall, 'b—^' ,SN,rmsbiasRt_sdma8, 'r- 

*' ,SN,rmsbiasRt_sdma4, 'g—o' ) 

xlabel (' Number of Arrays ',' Fontsize ', 12) 

ylabel( '\gamma,RMS(meters)' , 'Fontsize' ,12) 

legend (' LOS-only localization (fault angle) ', 'Proposed localization (1 
unknown reflector) ' , . . . 

'Proposed localization (2 unknown reflector) ') 

figure (7) 

plot(SN,rmsbiasRt_sdma5, 'b—^' , SN,rmsbiasRt_sdmalO, 'r-*') 
xlabel( 'Number Arrays' , 'Fontsize ' ,12) 
ylabel( '\gamma,RMS(meters)' , 'Fontsize' ,12) 

legend (' NLOS-only localization (1 reflector) ', 'NLOS-only localization 
(2 reflector)' ) 
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This is the Matlab code used in this thesis to evaluate the performance of the 
proposed source localization scheme in the presence one unknown reflector using the 
sequential least squares. 


%MODEL WHICH APPLIES THE SEQUENTIAL LEAST SQUARE IN THE LOCALIZATION 


clear; clc; clf reset; 

Var=0.25*pi''2/180''2; 

f=3E8; 

modulation 

TB=1/(4*f); 

= 4 

Nchips=64; 
seq 

SPC = 16; 
cp=3E8; 

lambda = cp/f; 
k = 2*pi/lambda; 

tend=SPC*Nchips; 

t=[0:(tend-1)]*TB/(tend-1); 

signal 

s=sin(2*pi*f*t); 

S=zeros(1,length(t)); 

S(l, :) =S(1, :)+exp(lj*(s)); 

%Signal Power 
P=10^-4; 

K=1.38*10^-23; 

Tk=400; 

B=3E4; 

SNR=10*logl0(P/(K*Tk*B)) 


% freq of baseband 
% time bandwidth product 
% Number of Chips / Time 
% Samples Per Chip 

% c/f 

% wavenumber 

% time snapshot end 
% discritization of 

%incident signal 


%% TARGET POSITION-REFLECTOR SIMULATION 
tic 

Xt=170; Yt=230; % targets 

position 

stdev=0.25; 

T=[Xt Yt]; 

Rt=sqrt(T(1)^2+T (2)^2); 
thol=10*pi/180; 

YRol=400; 

c=input( 'give the number of clusters-receivers:' ) 
array_size=input( 'give the size of the array:') 

N=input( 'give the number of sensors/cluster:') 
nruns=input( 'give the number of Monte carlo estimates:') 
it=input( 'give the number of iterations') 

Xs=zeros(l,c);Ys=zeros(l,c) ; 
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% Monte carlo simulation begin 


from 


for run=l:nruns; 
here 

fprintf( 'run:%2.Of\n' , run) ; 
YRol=YRol + stdev*randn; 

%Array positions 
for i=l:c; 

Ys(i)=45*rand; 

Xs(i)=45*rand; 

tan_thL=(Yt-Ys(i))/(Xt-Xs(i)) ; 
thL=atan(tan_thL); 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

000000000000000000000000000 

%Simulation of the reflector 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

000000000000000000000000000000 


XRsl= (1/ (1+tan (thol) '' 2 )) *Xs (i) + (tan (thol) / (1+tan (thol) ''2) ) * (Ys (i) - 
YRol); 

YRsl=tan(thol)*XRsl+YRol; 

Xsiml=2*XRsl-Xs(i) ; 

Ysiml=2*YRsl-Ys (i); 

thAl=2*thol+atan((Ysiml-T(2))/ (T (1)-Xsiml)); %DOA REF actual #1 
XR11=50; YRll=tan(thol)*XRll+YRol; 

XR12=200; YR12=tan(thol)*XR12+YRol; 
thll=atan((YRll-Ys(1))/(XRll-Xs(1)))+5*pi/180; 
thl2=atan((YR12-Ys(1))/(XR12-Xs(1)))-3*pi/180; 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

ooooooooooooooooooooooooooo 

%Signal strength simulation 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

ooooooooooooooooooooooooooooo 


%Array Gain 

GA=10*logl0(N)+2; %Assume dipole antenna elements 

%Free space path loss 

dL=sqrt((Xt-Xs(1))^2+(Yt-Ys(1))^2); 

LSL=-20*logl0(lambda)+20*1og10(dL)+21.98; 
dAl=sqrt((Xt-Xsiml)^2+(Yt-Ysiml)^2); 

LSAl=-20*logl0(lambda)+20*logl0(dAl)+21.98; 

%Reflection loss 

thTl=asin(sin(pi/2-thAl)/sqrt(5)) ; 

Href1=(cos(pi/2-thAl)-sqrt(5)*cos(thTl))/(cos(pi/2- 
thAl)+sqrt(5)*cos(thTl)); 

Lrefl=20*logl0(abs(Hrefl) ) ; 

%Received SNR 
SNRL=SNR-LSL+GA; 

SNRA1=SNR-LSA1+Lref1+GA; 


9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

ooooooooooooooooooooooooooooooooo 

%Received Signal- DOA Estimation 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

0000000000000000000000000000000000 

for j=l:it; 

[Bl,B2,B] =resignal(thL, thAl,N,array_size,k); 
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XL=B1*S; %Received signal 
XL=awgn(XL,SNRL, 'measured' ) ; 

XA1=B2*S; %Received signal 
XAl=awgn(XAl,SNRAl, 'measured' ); 

X=XA1+XL; 

thESTsdma=doa_sdma(X,B,N); 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

oooooooooooooooooooooooo 

%Assosiation 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

oooooooooooooooooooooooo 

REFl=find(thESTsdma<thll & thESTsdma>thl2); 

LOS=find(thESTsdma<thl2) ; 
thLest(i,j)=thESTsdma(LOS); 
thAest1(i,j)=thESTsdma(REFl) ; 
end 
end 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

0000000000000000000 

%Matrix Formulation for LOS and REF 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

0000000000000000000 

for j=l:it; 

for 1=1:c; 

Alos(i+c*j-c, :) = [! -tan(thLest (1,j))]; 

bios(i+c*j-c,1)=Ys(1)-tan(thLest(1,j))*Xs(1); 

Aref(i+c*j-c, :) = [! -tan(thAest1(1,j))]; 

bref (i+c*j-c,1)=Ys(1)-tan(thAest1(1,j))*Xs(1); 

end 

end 

Xestlos(:,1)=Alos(l:c,:)\blos(l:c,:); 

Xestref(:,1)=Aref(l:c,:)\bref(l:c,:); 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

0000000000000000000000000000000000000000000000000000 

%Reflector estimation using sequential least square 

9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 9 - 

000000000000000000000000000000000000000000000000000 

for j=l:c*(it-1) 

C (1: j+c-1,1:j+c-l) = (l/Var)*eye(j+c-1); 

Sigmalos(:,:)=inv(Alos(l:j+c-l,:)'*C(l:j+c-l,l:j+c-l)*Alos(l:j+c- 

1, : ) ) ; 

Klos(:,j)=(Sigmalos(:,:)*Alos(j+c,:)')/(Var+Alos(j+c,:)* 

Sigmalos(:, :)*Alos(j+c, :) ' ) ; 

Xestlos(:,j+1)=Xestlos(:,j)+Klos(:,j)*(blos(j+c,l)- 
Alos(j+c,:)*Xestlos(:,j)); 

Sigmaref(:, :)=inv(Aref(l:j+c-l, :) '*C(l:j+c-l,l:j+c-l)*Aref(1:j+c-1, : ) ) ; 

Kref ( : , j) = (Sigmaref ( : , : ) *Aref (j+c, -.)')/ (Var+Aref (j+c, : ) * 

Sigmaref(:, :)*Aref(j+c, :) ') ; 

Xestref(:,j+1)=Xestref(:,j)+Kref(:,j)*(bref(j+c,l)- 
Aref (j+c, :)*Xestref(:,j)); 

tan_thM=(Xestref(1, j)-Xestlos(1, j))/(Xestlos(2, j)-Xestref(2, j)) ; 
tho_est(j)=((pi/2)-atan(tan_thM)); 

YRo_est(j) = (Xestlos(l,j)+Xestref (1,j))/2- 

tan(tho_est(j))*((Xestlos(2, j)+Xestref(2, j))/2) ; 

for 1=1:c; 

XRsest ( j , 1) = (1 / (1+tan (tho_est (j) )''2) )*Xs(i) + (tan (tho_est ( j ) ) / (1+tan (tho 
_est (j) ) ''2) ) * (Ys (1) -YRo_est (j) ) ; 
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YRsest(j,i)=tan(tho_est(j))*XRsest(j)+YRo_est(j); 

Xsimest(j,i)=2*XRsest(j,i)-Xs(i); 

Ysimest(j,i)=2*YRsest(j , i)-Ys(i) ; 

Ar(i,:)=[l -tan(thLest(i))]; 

Ar(i+c, :) = [1 -tan(2*tho_est(j)-thAest1(i))] ; 
br (i,1)=Ys(i)-tan(thLest(i))*Xs(i) ; 

br(i+c,1)= Ysimest(j,i)-tan(2*tho_est(j)-thAest1(i))*Xsimest(j,i); 
end 

Xest2=Ar\br 

Rtest2=sqrt (Xest2 (1) ''2+Xest2 (2) ''2) ; 
biasRt_sdma2(run, j)=Rt-Rtest2 ; 
end 

%Los localization 
Xestl=Xestlos(:, 1) 

Rtestl=sqrt(Xestl(1)^2+Xestl (2)^2); 
biasRt_sdmal(run)=Rt-Rtestl; 
end 

meanbiasRt_sdmal=mean(biasRt_sdmal); 
position error 

varbiasRt_sdmal=var(biasRt_sdmal); 
of position error 

rmsbiasRt_sdmal=sqrt (meanbiasRt_sdmal. ''2+varbiasRt_sdmal) ; 
position error 

meanbiasRt_sdma2=mean(biasRt_sdma2); 
position error 

varbiasRt_sdma2=var(biasRt_sdma2); 
of position error 

rmsbiasRt_sdma2=sqrt(meanbiasRt_sdma2.^2+varbiasRt_sdma2); 
position error 

rmsbiasRt_sdmal=rmsbiasRt_sdmal*ones(1,8) 

SN=1:8; 
figure (1) 

plot(SN, rmsbiasRt_sdmal, 'b—^ , SN, rmsbiasRt_sdma2, 'r-*') 
xlabel (' Number of iterationsFontsize ', 12) 
ylabel( '\gamma,RMS(meters)' , 'Fontsize' ,12) 

legend (' LOS-only localization', 'LS proposed (known reflector) ') 


%Mean 

%Variance 

%RMS 


%Mean 

%Variance 

%RMS 
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This is the function used to determine the received signal in the array 

function [B1,B2,B] =resignal(thLOS,thREFl,N,array_size,k ) 


Ntheta=91; Nphi=3601; % Set up angle matrices 

theta = linspace(0,pi/2,Ntheta); 
phi = linspace(-pi,pi,Nphi); 

DTheta= theta(2)-theta(1); DPhi = phi(2)-phi(1); 

[Phi,Theta] = meshgrid(phi, theta) ; 

CT=cos(Theta) ; ST=sin(Theta) ; CP=cos(Phi); SP=sin(Phi); 

Rx=ST.*CP; Ry=ST.*SP; Rz=CT; % unit radial vector 

components 


min_spacing =array_size/50 ; 
sepn_min = 0; 

xp = array_size * rand(l,N); 
yp = array_size * rand(l,N); 


% in wavelengths (lambda) square grid 
% 1/lOth wavelength 
% initialize prior to 1st test 


Stats of Random Array 


while (sepn_min< min_spacing) % min spacing = .1 

wavelength 


rp = [xp;yp]'; 

sepn = pdist(rp); % dist between all 

elements in Random Array 

sepn_min = min(sepn); % find minimum separation 

zsepn = squareform(sepn); % convert linear sepn 

array to square 

[Rmin Cmin] = find(zsepn==sepn_min,1, 'first'); %returns index 
if (sepn_min< min_spacing) % then Move one of them 

xp(Rmin)=array_size * rand; 
yp(Rmin)=array_size * rand; 

end 


m=find(theta==pi/2) ; 
Rxm = Rx (m, : ) ; 

Rym = Ry (m, : ) ; 


% want row for xy plane 
% IxNphi vector 


% Calculate the actual received signal for each angle using the N 
element 

% receiver and sum all N signals, 
for n=l:N 

Bll(l,n) = exp(j*(k*(cos(thLOS)*xp(n)+sin(thLOS)*yp (n)))); 
B22(l,n) =exp(j*(k*( cos(thREFl)*xp(n)+sin(thREFl)*yp(n)))); 
end 


for n=l:3601; 

an(n)=(n/10-0.1)-180; 
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th(n)=an(n)*pi/180; 
for i=l:N 

B(i,n)=exp(lj*(k*(xp(i)* Rxm(n)+ yp(i)* Rym(n)))); 

end 

end 

B1=B11'; 

B2=B22 ' ; 


This is the function used to estimate the DOA of the incident signal using the SDMA 
receiver. 

function thESTsdma=doa_sdma(X,B, N) 

%% SDMA RECEIVER 
Nchips=64; 


SPC = 16; 

"5 

Samples Per Chip 

H = hadamard(Nchips); 

"6 

use Walsh Sequence 

codes=[]; 



for i=l:N 



codes(i, :)=H(i + 1, : ); 

"0 

don't want 1st row 

end 



% spread codes by SPC 



codes = repmat( codes,SPC,1); 

"5 

repeat by spread factor 

codes = reshape(codes,N,SPC*Nchips); 

"5 

move like chips together 

codes = codes*pi/2; 

"0 

pi/2 since inside exponential 


codes=codes'; 
r=codes*X; 

v=codes*B; 

R=conj(v)'*r; 

R=sum(R,2); 
maxR=max(abs(R)); 

Psdma=abs(R)/maxR; 

Psdma=[zeros(900,1) ; Psdma(902:2702, :); zeros(900, 1) ] ; 

%DOA SDMA 
Nphi=3601; 

phi = linspace(-pi,pi,Nphi); 
for i=l:Nphi 
if Psdma(i)>0.5; 

Psdma1(i)=Psdma(i) ; 

else 

Psdmal(i)=0 ; 

end 

end 

POSsdma = imregionalmax(Psdmal); % MATLAB rtn that finds 

local maxima->BW array "l"s 

PEAKsdma = find(POSsdma); % index pointers to local 

Maxima 
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[bs, is]=sort(Psdma(find(POSsdma) ),' descend ') ; % sort by highest peak- 
>iowest 

ThESTsdma= PEAKsdma(is); % index to peaks 

ThESTsdma= 180*phi(ThESTsdma)/pi; % Angle where peaks occur 

thESTsdma=ThESTsdma(1,1:2); 
thESTsdma=thESTsdma*pi/180; 


This is the function used to estimate the DOA of the incident signal using the MUSIC 
algorithm. 

doa_music.m 


function thESTmusic=doa_music(X,A,N) 

Nchips=64; % Number of Chips / Time seq 

SPC = 16; % Samples Per Chip 

tend=SPC*Nchips; 

Nphi=3601; 

phi = linspace(-pi,pi,Nphi); 

Rxx=X*X'/tend; 

[V, Dia]=eig(Rxx); 
matrix 

[Y,Index]=sort(diag(Dia)); 

EN=V(:,Index(1:N-2)); 
for i=l:3601; 

MU (i)=abs(A(:,i) '*EN*EN'*A(:,i)); 
end 

Pmu=l./MU; 
maxPmu=max(Pmu); 

Pmu=[Pmu(1802:3601) Pmu(1:1801)]; 

Pmusic=Pmu/maxPmu; 

Pmusic=[zeros(1, 900) Pmusic(901:2701) zeros(1, 900)] ; 

POSmusic = imregionalmax(Pmusic); 
local maxima->BW array "l"s 
PEAKmusic = find(POSmusic); 
local Maxima 

[bs,is]=sort(Pmusic(find(POSmusic) ),' descend' ); % sort by highest peak- 
>lowest 

ThESTmusic= PEAKmusic(is); 

ThESTmusic= 180*phi(ThESTmusic)/pi ; 
occur 

thESTmusic=ThESTmusic(1,1:3) ; 
thESTmusic=thESTmusic*pi/180; 


%autocorrelation matrix of signal 
%eigenvalues of autocorrelation 

%eigenvalues in descending order 
%noise subspace matrix 


%music pseudospectrum 


%MATLAB rtn that finds 
% index pointers to 


% index to peaks 
% Angle where peaks 
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